AD-A066  001  COAST  GUARD  WASHINGTON  0  C  OFFICE  OF  RESEARCH  AND  DC— ETC  F/6  17/7 

PRECISION  LORAN-C  NAVIGATION  FOR  THE  HARBOR  AND  HARBOR  ENTRANCE— ETC <U> 
HAY  BO  01  OLSEN*  J  M  LIOON*  A  J  SEOLOCK 
UNCLASSIFIED  USCS-D-3H-80  ML 


Report  No.  CQ-D-34-80 


irar 


PRECISION  LORAN-C  NAVIGATION  FOR  THE 


HARBOR  AND  HARBOR  ENTRANCE  AREA 


s 


DTIC 

|EI  P''TE1 
JUH2  7  1980 


MAY  1980 

Document  it  ivilliblt  to  tha  U.8.  public  through  tho 
National  Tachnical  Information  8arvlea. 
8prlngfiald.  Virginia  22161 


Prepared  by 

U.S.  DEPARTMENT  OF  TRANSPORTATION 
United  States  Coast  Guard 
Offica  of  Rotoarch  and  Davalopmant 
Washington,  D.C.  20593 


80  6  26  040 


4 


NOTICE 


This  document  is  disseminated  under  the  sponsorship  of  the  Department 
of  Transportation  in  the  interest  of  information  exchange.  The  United 
States  Government  assumes  no  liability  for  its  contents  or  use  thereof. 

The  contents  of  this  report  do  not  necessarilv  reflect  the  official  view 
or  policy  of  the  Coast  Guard;  and  they  do  not  constitute  a  standard, 
specification,  or  regulation. 

This  report,  or  portions  thereof  may  not  be  used  for  advertising  or 
sales  promotion  purposes.  Citation  of  trade  names  and  manufacturers 
does  not  constitute  endorsement  or  approval  of  such  products. 


2.  Government  Accession  No. 

P&'hCPrL  <t><t>±- 


1  »"■ . .  HI|MH 

Depart ment-of  Transportation 
tfrS.  Coast  Guard 

Office  of  Research  and  Development  (G-BST-1) 
Washington,  DC  20593 


Technical  Report  Documentation  Page 


3.  Recipient  s  Catalog  No. 


May  08 


6.  Performing  Orgoni  option  Code 


I.  Performing  Organisation  Report  No. 

CG-D-34-8(/ 


JO  Work  Unit  HU.  1TRAIS) 


13.  Typ*  ol  Rtp.fl  and  Period  Ctvtrtl 


Id.  Saantonng  Aptncy  Cud* 

G-DST-1 


12.  Spofitoripp  Aptney  Nam.  and  Addratt 

Department  of  Transportation 
U.S.  Coast  Guard 

Office  of  Research  and  Development  (G-DST-1) 
Washington,  DC  20593 


1$  Supplementary  Note* 

This  report  updates  and  consolidates  three  papers  presented  at  the  Eighth  Annual  Technical 
Symposium  of  the  Wild  Goose  Association  in  Williamsburg,  Virginia,  October  17-19,  1979.  The 
papers  will  be  published  as  part  of  the  proceedings  of  that  symposium. 


16.  Abstract 

'  The  U.  S.  Coast  Guard  is  engaged  in  a  long-term  program  to  evaluate  the  potential  and  develop  the  technology  for 
precision  navigation  using  Loran-C.  The  Coast  Guard  has  taken  a  three-pronged  approach  to  harbor  and  harbor  entrance 
(HHEl  navigation— user  equipment  development,  time  difference  (TD)  surveying,  and  grid  stability  analysis.  These  three 
aspects  of  HHE  Loran-C  navigation  were  addressed  in  separate  papers  presented  in  October,  1979.  This  report  updates  and 
consolidates  the  three  papers,  and  is  intended  to  serve  as  a  source  document  on  the  Coast  Guardis  approach  to  HHE 
Loran-C  navigation. 

User  equipment  is  being  developed  which  combines  time  differences  from  a  precision  Loran-C  receiver  with  survey 
information  from  a  data  tape  cassette  using  simple  algorithms  to  present  the  piloting  situation  on  a  CRT  display.  The 
algorithms,  hardware,  and  software  are  described,  and  the  results  of  summer  and  winter  1979  tests  on  the  St.  Marys  River 
are  given.  It  is  concluded  that  accuracies  of  better  than  25  meters  are  achievable. 

Three  methods  are  described  for  surveying  the  Loran-C  time  difference  coordinates  of  waypoints  and  prominent 
physical  features  of  a  harbor  area  for  precisian  Loran-C  navigation.  Results  are  presented  from  field  measurements  on  the 
St.  Marys  River  Mini  Loran-C  Chain  using  existing  visual  aids  to  navigation  as  a  position  reference.  The  survey  methods 
described  provide  for  accurate  and  timely  Loran-C  grid  survey  over  a  wide  range  of  harbor  areas.  The  general  features  of 
each  method  and  details  of  the  Visual  aids  to  navigation  approach  are  presented.  Analysis  of  field  measurements  show 
survey  capability  on  the  order  of  ten  nanoseconds.  . 

Temporal  instabilities  on  the  order  of  several  hundred  nanoseconds  were  measured  over  the  service  area  of  the  St. 
Marys  River  Mini-Chain  between  Pall  1977  and  Spring  1978.  These  instabilities,  if  true,  could  render  the  mini-chain 
unusable  for  precision  navigation  on  the  St.  Marys  River.  A  careful  recollection  effort  was  deemed  necessary  to  verify  the 
chain's  performance  and  to  definitively  judge  Its  navigation  capabilities.  In  addition,  a  methodology  is  desired  for 
establishing  the  stability  of  Loran-C  signals  in  the  general  harbor  and  harbor  entrance  environment.  This  methodology  is 
being  investigated  using  the  mini-chain.  Time  difference  data  is  being  collected  at  three  fixed  monitor  sites,  and  is  beiiy 
analysed  by  a  rather  simple  model  which  separates  variations  into  uniform  propagation  effects,  local  errors  and  control 
errors.  The  propagation  component  extracted  from  the  data  exhibits  a  high  correlation  with  temperature,  which  can  be 
attributed  primarily  to  the  vertical  lapse  rate  of  the  index  of  refraction. 
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PREFACE 


The  U.  S.  Coast  Guard  is  engaged  in  a  long-term  program  to  evaluate  the  potential 
and  develop  the  technology  for  precision  navigation  using  Loran-C.  The  basic  method 
is  to  use  Loran-C  in  the  repeatable  mode  for  real  time  2  or  3  line-of-position  precision 
navigation,  or  piloting,  in  the  vicinity  of  surveyed  positions.  This  method  can  be 
augmented  by  offset  or  proportional  differential  corrections.  The  Coast  Guard  has 
taken  a  three-pronged  approach  to  harbor  and  harbor  entrance  (HHE)  navigation— user 
equipment  development,  time  difference  (TD)  surveying,  and  grid  stability  analysis. 
These  three  aspects  of  HHE  Loran-C  navigation  were  addressed  in  separate  papers - 
presented  in  October,  1979.  This  report  updates  and  consolidates  the  three  papers,  and 
is  intended  to  serve  as  a  source  document  on  the  Coast  Guard's  approach  to  HHE 
Loran-C  navigation. 

Section  1  presents  the  paper  by  Ligon  and  Edwards*  on  user  equipment  develop¬ 
ment.  The  PILOT  (Precision  Intracoastal  Loran  Translocator)  equipment  described 
represents  an  effort  to  produce  a  com m  er cially-reali  zable  precision  Loran-C 
navigator.  A  basic  philosophy  in  PILOT'S  design  is  the  minimization  of  real-time 
computations.  Many  of  the  required  parameters  are  pre -computed  and  stored  on  a 
cassette  tape  which  includes  as  well  the  definition  of  the  channel,  navaids  and 
geographic  features  of  a  waterway.  Included  in  this  section  are  the  results  of  a  brief 
evaluation  of  PILOT  conducted  on  the  St.  Marys  River  during  October,  1979.  Several 
additional  sets  of  this  equipment  are  presently  being  produced  for  a  1980-1981 
evaluation  by  various  government  and  commercial  interests.  Appendices  A-I  present 
the  development  of  the  algorithms  used  within  PILOT. 

Section  2  presents  the  papa1  by  Sedlock^  on  time  difference  surveying.  A  survey 
is  required  to  determine  the  repeatable  Loran-C  TDs  which  define  a  waterway,  and 
which  are  encoded  onto  cassette  tapes  for  PILOT.  A  Time  Difference  Survey  Set 
(TDSS)  and  accompanying  statistical  algorithms  provide  the  means  for  conducting  the 
survey.  Several  examples  show  survey  data  gathered  on  the  St.  Marys  River. 

Section  3  concludes  the  report  with  a  paper  on  grid  stability  analysis  by  Olsen  and 
Isgett3.  Grid  stability  is  fundamental  to  precise  navigation  using  repeatable  Loran-C 
coordinates.  The  paper  concentrates  on  the  St.  Marys  River  Mini-Loran-C  Chain,  and 
mentions  the  extension  of  the  instrumentation  and  analytical  techniques  to  other  HHE 
areas  around  the  country. 
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Section  1 


USER  EQUIPMENT  DEVELOPMENT 

LORAN-C  PILOTING 
Introduction 

As  a  navigator  proceeds  from  the  Coastal  Confluence  Zone  into  harbor  entrances, 
harbors  and  rivers,  his  requirements  for  both  positioning  accuracy  and  timeliness 
increase  on  the  order  of  500  times.  The  number  of  information  sources  also  increases 
and  the  navigator  must  integrate  the  data  from  all  of  these  sources.  He  has  crossed 
the  boundary  separating  piloting  from  navigation. 

Large  areas  of  the  U.S.  Coastal  Confluence  Zone  (CCZ)  are  presently,  or  will  shortly 
be,  covered  by  Loran-C  signals  which  are  charted  to  an  accuracy  of  0.4  km.  It  is 
known  that  these  signals  are  repeatable  to  within  a  few  tens  of  meters.  These  signals 
are  also  present  in  harbors  and  harbor  entrances  and  potentially  form  the  basis  for  a 
precision  navigation  system.  This  section  discusses  a  method  for  exploiting  these 
Loran-C  signals  and  describes  the  PILOTing  hardware  which  is  the  keystone  for  this 
exploitation.  The  name  PILOT— Precision  Intracoastal  Loran  Translocator— is  not 
meant  to  imply  automatic  piloting,  but  rather  to  highlight  the  significantly  different 
frame  of  reference  which  exists  in  the  piloting  domain. 

Loran-C  PILOTing 

The  first  two  words  forming  the  acronym  "PILOT"  refer  to  the  environment;  the  last 
two  refer  to  the  method. 

Intracoastal  waters  are  typically  characterized  by  visible  or  physical  features:  channel 
boundaries,  visual  aids-to-navigation,  coastlines.  They  are  also  characterized  by 
significant  overland  Loran-C  propagation  paths  and  warped  time  difference  grids.  In 
these  waters,  positioning  in  an  absolute,  geodetic  sense  (latitude/longitude)  is  of  little 
value.  Also,  the  orientation  of  the  vessel  is  as  important  as  its  position.  Precision 
must  be  relative  to  these  intracoastal  realities— "precise"  enough  for  safe,  speedy 
vessel  operation. 

The  method  is  basically  Loran-C  translocation.  Given  a  reference  site  and  its 
corresponding  Loran-C  time  differences,  estimate  the  location  of  a  Loran-C  receiver 
at  another  site  in  the  vicinity  of  the  reference  site,  based  upon  the  time  difference 
measurements  taken  there.  This  is  essentially  the  same  as  the  method  used  in 
conventional  satellite  position  translocation  and  in  many  loran  waypoint  navigation 
systems. 
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ENVIRONMENTAL  REQUIREMENTS 
Calibration 

If  the  latitude  and  longitude  of  all  the  essential  navigation  features  in  a  particular 
river  or  harbor  were  known,  and  if  there  were  a  perfect  propagation  prediction  model, 
then  calibration  would  be  unnecessary.  There  are  over  100,000  visual  navigation  aids 
installed  in  the  U.S.  intracoastal  waters,  and  the  latitude  and  longitude  are  known  for 
very  few  of  them;  yet,  they  serve  their  purpose.  In  the  past,  attempts  have  been  made 
to  utilize  the  most  accurate  Loran-C  prediction  methods  in  the  piloting  environment, 
but  unfortunately  these  require  a  larger  geological  data  base  than  the  corresponding 
loran  data  base  for  river  and  harbor  "calibration"  by  measurement,  and  still  do  not 
provide  sufficient  accuracy.4 

When  calibration  is  by  measurement,  coding  delays  or  System  Area  Monitor  reference 
time  differences  are  held  constant  and  the  time  difference  field  is  measured  as  it  is. 
No  attempt  is  made  to  adjust  the  coding  delays  to  cover  the  service  area  in  some  sort 
of  "least  squares"  error  sense  to  fit  the  imperfect  propagation  model.  Calibration  by 
measurement  brings  up  five  questions:  the  questions  of  density  of  survey,  of 
coordinate  conversion  algorithm,  of  positioning,  of  survey  time  difference 
measurement,  and  of  temporal  stability.  The  questions  of  density  of  survey  and 
coordinate  conversion  algorithm  are  very  tightly  coupled.  On  the  St.  Marys  River,  a 
survey  interval  of  approximately  3  kilometers  coupled  with  navigational  aid  relative 
position  taken  from  the  river  charts  and  processed  using  the  PILOT  coordinate 
conversion  algorithm  provided  "sufficient"  navigational  accuracy.  That  is,  the  PILOT 
showed  the  vessel  to  be  in  the  center  of  the  channel  and  passing  the  proper 
navigational  aids  when  it  appeared  to  the  observer  that  this:  was  the  case.  The 
remaining  three  questions  have  been  addressed  and  are  discussed  in  greater  detail  in 
Sections  2  and  3,  and  in  the  paper  by  Johler5.  In  summary,  there  is  a  very  effective 
method  for  surveying  the  Loran-C  time  differences  relative  to  the  visual  aids  and 
features  used  by  the  pilot  himself,  and  the  stability  question  is  vital  and  is  ultimately 
solved  by  a  form  of  differential  corrections. 

Reliability  and  Speed 

The  piloting  environment  places  constraints  upon  navigational  equipment  that  do  not 
exist  in  the  more  open  waters.  The  pilot  needs  much  higher  accuracy  and  reliability. 
The  safety  of  his  vessel  cannot  tolerate  the  loss  of  a  fix  caused  by  failure  of  one  of  the 
navigational  transmitting  stations.  In  Loran-C  this  translates  to  3-LOP  fixes,  with 
drop-back  2-LOP  fixes,  including  those  when  the  Master  station  has  a  casualty. 

The  pilot  does  not  have  time  to  convert  numbers  to  a  position  on  a  chart.  His 
navigational  display  should  provide  him  with  an  immediate  picture  of  his  situation  so 
that  he  can  rapidly  respond  to  it.  This  translates  to  a  visual  plan  display  with 
navigational  aids,  channel  boundaries  and  some  shoreline  features  indicated. 

Other  Issues 

It  is  possible,  and  PILOT  is  the  demonstration  case,  to  meet  the  requirements  of  the 
piloting  environment  with  a  compact  user  equipment  for  approximately  $20,000, 
including  the  cost  of  a  "better  grade"  Loran-C  receiver.  The  pilot  will  probably  not  be 
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able  to  obtain  satisfactory  operations  using  a  low-cost  CCZ  Loran-C  receiver. 
PILOT 


Description 

The  PILOT,  developed  for  the  U.S.  Coast  Guard  by  the  Applied  Physics  Laboratory  of 
Johns  Hopkins  University,  is  an  electronic  aid  for  piloting  vessels  in  harbors  and 
rivers.  Prerecorded  tape  cartridges,  containing  a  sequence  of  chartlets  and  other 
navigation  information,  provide  the  PILOT  with  a  degree  of  "local  knowledge."  The 
vessel's  present  position  and  heading,  continuously  determined  from  a  loran  receiver 
and  the  vessel's  gyro,  are  displayed  on  the  current  area  chartlets.  Chartlets  of  two 
different  scales  (master  and  detail)  are  always  available  for  operator  selection. 
Position  information  relative  to  waypoints  (intermediate  destinations)  is  displayed  to 
the  left  of  the  chartlet.  A  horizontal  bar  graph,  representing  the  vessel's  relative 
cross  track  position,  can  be  displayed  along  the  bottom  of  the  CRT.  A  predictor  option 
displays  the  ship's  future  position,  based  upon  present  course  made  good  and  turning 
rate,  for  any  operator  selectable  time  period. 

There  have  been  many  systems  constructed  for  similar  purposes  in  the  past.  The 
earliest  was  probably  the  Decca  track  plotter  in  the  early  1960s.  Modem  track  plotters 
are  found  on  many  fishing  vessels.  Table  1-1  lists  some  of  these  coordinate 

TABLE  1-1.  COMPARISON  OF  VARIOUS  TYPES  OF  USER  EQUIPMENT 
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X-Y 
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mm 
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Video 
Graphics 
t  Alpha- 
numerics 

Shows  vessel  position  and 
orientation.  Ultimately 
required  Intensive  survey. 
Primarily  an  R/D  device. 
Very  long  chartlet  change 
time. 

Modified  User  11 

Inverse 

Prediction 

THt 

Color 
Graphics 
&  Alpha- 
numerics 

As  desired. 
Data  on  disc 
memory. 

Color  graphics  and 

RAYDIST  and  Miniranger 
sensors  added. 

TOT  " 

6-Matrix 
with  flat 
earth  hyper¬ 
bolic  grid 

i 

Video 

graphics 
chartlets 
analog,  A 
alpha- 
nuaerlcs 

As  desired. 

Display  and 

conversion 
data  on  tape 
cassette. 

my 

conversion/display  devices  and  comments  upon  their  features  as  they  relate  to  the 
harbor  and  river  environment.  As  can  be  seen,  the  PILOT  is  a  logical  progression  from 
these  devices.  It  solves  the  coordinate  conversion  problem  which  had  in  the  past  been 
the  major  obstacle  to  HHE  implementation,  and  it  demonstrates  the  feasibility  of 
reliable  PILOTing  using  a  single,  compact  video  display  terminal. 

Navigational  Charts 

The  PILOT  is  basically  a  data  processor.  All  coordinate  conversion  constants  and 
navigational  background  information  are  provided  on  a  cassette  which  is,  in  effect,  a 
Loran-C  navigational  chart  on  a  magnetic  tape  medium.  A  single  cassette  can  contain 
as  many  as  150  to  200  of  these  chartlets— more  than  enough  to  permit  a  complete 
transit  of  the  St.  Marys  River,  which  is  approximately  100  kilometers  long.  Section  2 
contains  a  deeper  discussion  of  the  chartlet  calibration.  In  summary,  the  method 
involves  the  creation  of  a  linked  mosaic  of  ideal  flat-earth  hyperbolic  Loran-C  grids 
fixed  at  the  calibration  points,  with  aids  and  geographic  features  charted  onto  these 
grids.  This  implies  a  slight  distortion  of  the  physical  world  if  the  real  Loran-C  grid  is 
distorted.  On  the  St.  Marys  River,  this  distortion  was  small  enough  to  be  unnoticeable, 
even  with  calibration  point  separations  as  great  as  10  kilometers.  Had  the  distortion 
become  too  great,  additional  survey  points  would  have  been  necessary.  Since  the 
identical  model  is  used  within  the  PILOT  for  coordinate  conversion,  the  charting 
accuracy  is  as  good  as  the  repeatable  accuracy. 

HARDWARE  DESCRIPTION 

Graphics  Terminal 

The  nucleus  of  the  PILOT  system  (Figure  1-1)  is  a  Hewlett  Packard  2649A 
microprogrammable  graphics  terminal.  This  OEM  device  was  selected  because  it  has  a 
separate  graphics  processor  with  memory,  dual  tape  cartridge  units,  and  an  8080 
microprocessor  that  could  be  modified  and  programmed  as  required.  Modifications  to 
the  HP-2649A  terminal  included:  converting  the  8080  microprocessor  from  software 
math  to  hardware  math  by  adding  an  Advanced  Micro  Devices  (AMD)  9511  arithmetic 
processing  unit,  developing  a  two-receiver  interface  board  to  mount  inside  the 
terminal,  developing  an  interface  board  to  connect  to  the  ship's  gyro  and  to  a  time 
difference  bias  (differential  correction)  box,  replacing  the  large  general  purpose 
keyboard  with  a  small  predefined  keypad,  and  building  a  short  base  for  the  terminal  to 
serve  as  a  cable  junction  box.  The  system  block  diagram  is  shown  in  Figure  1-2. 

Receiver  Interface  Board 

The  dual  receiver  interface  board  is  designed  around  an  8085  microprocessor  and  a 
second  AMD  9511  math  unit.  This  microprocessor  performs  the  following  functions: 
input,  convert,  identify,  edit,  filter  and  dead  reckon  (DR).  The  use  of  two  loran 
receivers  permits  cross  chain  position  fixes.  The  data  format  of  most  receivers  can  be 
accommodated  (up  to  50  parallel  lines  per  receiver)  by  changing  the  input  software 
module.  A  serial  receiver  interface  board  has  been  designed,  and  will  be  used  to 
operate  PILOT  with  the  forthcoming  Internav  model  LC-404  Loran-C  Navigation  and 
Monitor  Receiver,  a  microprocessor  version  of  the  presently-used  Internav  MK  in 
receiver. 
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FIGURE  1-1.  PILOT 


FIGURE  1-2.  PILOT  SYSTEM  BLOCK  DIAGRAM 
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FIGURE  1-3.  PILOT  FUNCTIONAL  DATA  FLOW  DIAGRAM 
Heading  and  TD  Bias  Board 

The  heading  and  TD  bias  board  includes  a  10-bit  synchro-to-digital  converter  and  can 
accept  heading  information  from  either  a  gyrocompass  or  a  magnetic  compass  with  a 
shaft  encoder.  Storage  registers  on  this  board  can  store  three  sets  of  TD  bias  numbers 
(+  7999  nanoseconds).  These  registers  can  be  loaded  either  by  the  operator  from  the 
keypad  or  from  a  remote  box  or  modem  via  an  interface  cable. 

Keypad 

Operation  of  the  PILOT  terminal  was  simplified  by  replacing  the  HP  keyboard  (105 
keys)  with  a  small  keypad  (23  keys)  having  eleven  function  keys  plus  twelve  numerical 
and  cursor  control  keys. 

SOFTWARE  DESCRIPTION 

Figure  1-3  is  the  PILOT  functional  data  flow  diagram. 

Memory 

Approximately  17,000  bytes  of  machine  language  code  were  developed  at  APL  for  the 
PILOT  terminal,  and  an  additional  40,000  bytes  of  the  original  HP  code  was  retained. 
Structured  programming  and  assembly  language  was  used  for  maximum  efficiency. 

TD  Filtering 

Digital  filtering  is  used  on  each  TD  as  part  of  the  preprocessing  performed  on  the 
receiver  interface  board.  An  alpha/beta  filter  was  modified  to  induct  turning 
acceleration  feedback.  For  receivers  with  relatively  slow  data  rates  (i.e.,  less  than 
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two  sets  of  TDs  per  second),  the  TDs  are  dead-reckoned  every  100  milliseconds  using 
the  last  known  velocity  and  turning  acceleration.  Three  different  filter  time  constants 
may  be  selected  by  the  operator  from  the  keypad.  The  development  of  PILOTS  filter 
algorithms  is  given  in  Appendices  D-H. 

Coordinate  Conversion 

The  8080  microprocessor  performs  a  full  order,  iterative  transformation  approximately 
twice  per  second.  When  three  good  TDs  are  available,  a  2x3  minimum  variance  "G" 
matrix  is  used.  Each  term  of  the  matrix  is  preweighted  to  produce  a  best  fit  with  the 
expected  relative  signal  strength.  When  only  two  TDs  are  available,  one  of  three  2x2 
unweighted  matrices  is  used.  The  development  of  the  "G"  matrix  and  the  iterative 
exact  solution  to  the  non-linear  loran  equations  are  presented  in  Appendices  A-C. 

Chartlet  Cassettes 

Chartlet  cassettes  contain  an  index  file,  master  files  and  detail  files.  The  index  file 
contains  a  title  block  and  a  list  of  all  master  ehartlets  on  the  cassette.  Each  master 
file  contains  the  graphics  for  8-16  miles  of  track,  the  area  matrix  coefficients, 
transmitter  coordinates  and  supplemental  data  such  as  display  origin,  scale,  rotation, 
and  geometric  dilution  of  precision  (GDOP— developed  in  Appendix  I).  Each  detail  file 
contains  the  graphics  for  1-2  miles  of  track,  the  TDs  and  x-y  coordinates  of  the 
current  waypoint,  bearing  angles  to  and  from  the  WP  and  supplemental  data. 

Master  ehartlets  provide  "look  ahead"  by  showing  the  next  several  waypoints;  detail 
ehartlets  provide  a  closer  view  of  the  vessel's  current  situation.  Each  master  file  is 
followed  by  one  or  more  detail  files.  North-up  or  track-ip  ehartlets  may  be  used. 

Chartlets  can  be  developed  using  a  digitizer,  minicomputer  or  calculator,  plus  the 
PILOT  terminal  reconfigured  as  a  "stock"  HP-2649A  terminal.  The  original  chartlet 
cassettes  for  the  St.  Marys  River  were  produced  at  APL,  and  the  Coast  Guard 
Research  and  Development  Center  is  currently  setting  up  a  system  for  producing  them 
for  other  locations. 

INITIAL  TEST  RESULTS 

The  PILOT  system  was  initially  tested  on  the  St.  Marys  River  during  the  first  week  in 
October  1979. 

The  Receiver 

The  Loran-C  receiver  used  was  an  Internav  MK  HI  Survey  Receiver.  The  receiver  had 
been  previously  tested  in  the  laboratory  and  could  typically  be  expected  to  provide  an 
accuracy  of  30  nanoseconds  under  a  variety  of  signal  conditions.  The  St.  Marys  River 
signal  environment  was  so  benign  that  performance  in  the  10  nanosecond  range 
(relative  to  a  reference  Austron-5000  monitor  receiver)  was  observed. 

The  PILOT 

Overall  performance  was  excellent.  The  PILOT  successfully  combined  the  receiver  TDs 
with  survey  data  taken  during  August  and  September  and  chartlets  prepared  from  the 


NOAA  navigational  charts  to  provide  accurate  position  and  situation  information.  The 
references  in  most  areas  were  visual  ranges  and  buoys,  and  were  quite  good.  The 
navigational  solutions  were  based  upon  waypoint  surveys,  and  the  chartlets  were 
calibrated  by  using  the  surveyed  location  of  a  single  aid  on  each  one.  Quantifiable 
accuracy  tests  will  be  completed  next  year,  but  visual  accuracy  was  excellent.  There 
was  no  noticeable  grid  warp— possibly  the  result  of  the  low  conductivity  of  both  the 
fresh  water  and  the  Laurentian  Shield  upon  which  the  area  is  located.  A  repeatable 
mode  error  budget  is  given  in  Section  3. 

WINTER  DEPLOYMENT 

Following  the  preliminary  field  tests,  the  PILOT  system  was  installed  aboard  the  U.S. 
Coast  Guard  Cutter  Katmai  Bay  for  operational  deployment  during  the  winter  of 

1979- 1980.  Receiver  and  PILOT  operation  was  easily  mastered  by  all  members  of  the 
ship's  bridge  crew.  The  display  of  alphanum  erics,  graphics  and  the  cross-track  distance 
bar  offered  something  for  each  of  the  various  levels  of  interest,  from  commanding 
officer  to  helmsman.  PILOT  was  shown  to  be  useful  in  thick  fog  and  for  starting  turns 
based  upon  the  distance  to  go.  Though  the  expiration  of  the  winter  navigation  program 
limited  icebreaking  activities,  the  consensus  is  that  PILOT  should  be  a  useful  aid  in 
breaking  channel  edges  not  marked  with  buoys. 

THE  FUTURE 

A  data  logging  capability  is  planned  f or  PILOT,  whereby  a  record  of  a  voyage  will  be 
recorded  on  the  second  tape  drive  of  the  HP-2649A  terminal.  The  remaining 
development  tasks  include  interfacing  to  the  LC-404  Loran-C  receiver  and  installing 
the  algorithms  to  allow  a  fall-back  Master  independent  2-LOP  capability. 

PILOT  will  be  reinstalled  aboard  the  Katmai  Bay  in  June  1980  for  further  evaluation 
and  for  demonstration  to  commercial  shipping  interests.  Several  PILOT  systems  have 
been  offered  to  the  carriers  for  evaluation  aboard  their  vessels  in  the  St.  Marys  River 
and  cm  the  Great  Lakes  between  August  and  December  1980.  These  deployments  will 
be  followed  by  quantitative  and  human  factor  tests. 

A  "stripped-down"  version  of  PILOT,  designated  PLAD  for  Portable  Loran  Assist 
Device,  is  presently  under  development  at  the  Applied  Physics  Laboratory.  PLAD  will 
use  a  simple  digital  display  in  a  hand-held  data  terminal  in  lieu  of  the  graphical  display 
of  the  HP-2649A,  and  will  store  waypoint  information  in  semiconductor  memory.  The 
device  is  planned  for  initial  testing  on  the  Delaware  Bay  during  the  winter  of 

1980- 1981. 

CONCLUSION 

The  PILOT  opens  an  entirely  new  area  for  precision  Loran-C  navigation. 
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Section  2 


TIME  DIFFERENCE  SURVEYING 


INTRODUCTION 

The  Coast  Guard  is  currently  engaged  in  a  development  program  for  precision  Loran-C 
navigation  in  harbor  areas.  The  key  element  in  this  program  is  the  PILOT  user 
equipment  described  in  Section  1.  PILOT  provides  position  estimates  based  upon 
Loran-C  time  difference  (TD)  measurements,  a  flat-earth  hyperbolic  grid  (FEHG) 
coordinate  conversion  algorithm,  and  ehartlets  stored  on  a  magnetic  tape  cartridge. 
The  ehartlets  contain  position  and  TD  coordinates  for  waypoints,  a  description  of 
channel  boundaries  and  aids  to  navigation,  and  computational  constants  for  the  FEHG 
algorithm.  If  harbor  areas  were  accurately  surveyed  (in  the  geodetic  sense)  and 
Loran-C  propagation  models  were  highly  accurate,  the  data  necessary  to  produce 
ehartlets  could  be  directly  calculated.  Since  neither  is  true,  a  technique  has  been 
developed  to  survey  the  TD  coordinates  for  waypoints  and  to  test  for  TD  grid  warp 
which  could  affect  the  accuracy  of  the  navigation  solution.  This  survey  technique 
applies  to  a  wide  range  of  navigation  scenarios  including  relatively  wide  open  areas 
such  as  San  Francisco  Bay  and  severely  restricted  channels  such  as  the  St.  Marys 
River.  In  addition,  the  survey  technique  is  cost  effective,  and  data  reduction  and 
analysis  can  be  accomplished  in  the  field  as  the  survey  progresses. 

CLASSIC  APPROACH 

The  classic  approach  to  the  survey  problem  is  to  use  a  high  accuracy  reference  system 
to  position  a  survey  vessel  at  a  point  and  simultaneously  record  position  and  TD 
coordinates.  This  is  the  basic  approach  that  was  used  in  initial  attempts  to  survey 
waypoints  on  the  St.  Marys  River.  Several  difficulties  arise  with  this  approach.  An 
accurate  geodetic  description  of  channel  boundaries  and  aids  to  navigation  does  not 
exist,  and  where  accurate  coordinates  do  exist,  they  are  not  always  consistent  with 
geodetic  control  points  ashore.  Often  existing  geodetic  control  points  ashore  are  not 
recoverable  or  are  in  locations  unsuitable  for  locating  the  reference  system 
transponders,  and  additional  control  points  must  be  surveyed.  The  cost  of  establishing 
the  necessary  geodetic  control  ashore,  operating  the  position  reference  system,  and 
performing  the  TD  survey  are  substantial.  As  the  navigation  channel  becomes  more 
restrictive  the  inconsistency  between  survey  coordinates  and  coordinates  for  channel 
boundaries  and  aids  to  navigation  becomes  intolerable.  The  costs  both  in  time  and 
dollars  for  a  geodetic  survey  of  channel  boundaries  and  aids  to  navigation  for  an  entire 
harbor  are  totally  prohibitive. 

A  NEW  APPROACH 

The  basic  problem  producing  ehartlets  for  precision  Loran-C  navigation  is  relating  the 
Loran-C  TD  grid  to  the  navigation  situation  of  channel  boundaries,  shoals, 
aids-to-navigation,  etc.  This  section  describes  an  approach  which  has  been  developed 
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for  determining  TD  coordinates  of  the  physical  features  of  a  harbor  by  surveying.  In 
addition,  analysis  tools  are  used  to  obtain  position  coordinates  which  describe  the 
spatial  orientation  of  the  physical  features  consistent  with  the  TD  grid.  Although 
these  position  coordinates  are  accurate  enough  in  a  local  sense  that  deviations  from 
geodetic  coordinates  are  imperceptible,  no  attempt  is  made  to  tie  the  position 
coordinates  to  a  geodetic  reference. 

Since  visual  aids-to-navigation  form  the  basic  position  reference,  the  approach  has 
been  termed  the  Visual  Grid  Survey.  A  TD  measurement  set  has  been  developed  by  the 
Coast  Guard  Research  and  Development  Center  to  accurately  and  efficiently  measure, 
record  and  process  TD  information.  Analysis  tools  include  the  techniques  for 
determining  waypoint  TDs  and  x-y  coordinates  using  the  FEHG  algorithm  developed  for 
the  PILOT  user  equipment.  The  FEHG  algorithm  calculates  position  relative  to  a 
reference  waypoint  based  upon  the  difference  between  observed  TDs  and  the  TDs  for 
the  reference  point. 

General  Procedure 

The  general  procedures  for  a  Visual  Grid  Survey  are; 

a.  Choose  waypoints  from  harbor  charts  and  estimate  the  waypoint  coordinates 
(e.g.  latitude,  longitude). 

b.  Compute  the  FEHG  parameters  for  these  positions. 

c.  Estimate  the  positions  of  surrounding  features  needed  for  PILOT  chartlets 
relative  to  the  waypoint,  i.e.  construct  first  cut  chartlets. 

d.  Survey  the  TDs  of  ranges,  shoals,  channel  edges,  aids-to-navigation,  etc.  in  the 
area  of  each  waypoint. 

e.  Define  the  TDs  of  each  waypoint. 

f.  Refine  the  actual  position  offsets  of  charted  features  with  respect  to  each 
waypoint  using  the  FEHG  algorithm  and  channel  edge  TD  data. 

g.  Link  the  chartlets  from  a  central  or  major  waypoint  by  calculating  the  position 
of  adjacent  waypoints  from  the  surveyed  waypoint  TDs  and  the  FEHG  algorithm. 

This  procedure  results  in  chartlets  that  are  locally  exact  models  of  reality.  Charting 
and  chart-piokoff  errors  are  removed  but  local  grid  warp  and  TD  survey  errors  remain. 

Waypoints 

As  outlined  above,  the  first  step  in  the  survey  chartlet  preparation  process  is  the 
selection  of  waypoints.  The  waypoints  are  defined  based  upon  navigation  charts  and 
the  knowledge  of  traffic  patterns.  The  intersection  of  two  ranges,  intersection  of  two 
channel  centerlines,  intersection  of  right  half  channel  centerlines,  or  the  intersection 
of  two  commonly  used  tracklines  all  define  possible  waypoints.  There  are  two  general 
approaches  to  determining  waypoint  TDs.  The  first  approach  is  used  where  the 
waypoint  is  defined  by  the  intersection  of  two  visual  ranges.  The  visual  ranges  provide 


very  precise  crosstrack  information.  In  a  small  region  near  the  waypoint  the  Loran-C 
TD  grid  can  be  approximated  by  the  linear  model  below: 


TDX-TDX0  =  an(x-x0)  +  ai2(y-y0)  +  nx  (2-1) 

TDY-TDY0  =  a2i(x-xQ)  +  a22(y-yQ>  +  ny  (2-2) 

TDZ-TDZ0  =  a3i(x-Xo)  +  a32(y-yo)  +  nz  (2-3) 


where 


TDX,  TDY,  TDZ  are  observed  TDs 

TDX0,  TDY0,  TDZ0  are  waypoint  TDs 

x,  y  are  position  coordinates  of  the  observation 

xQ,  y0  are  waypoint  position  coordinates  of  the  waypoint 

ajj  are  coefficients  of  the  gradient  matrix  (directional  derivatives  of  the  TD  grid) 
nx,  ny,  nz  are  error  terms  due  to  noise  and  nonlinearity 

Along  the  centerline  of  a  visual  range, 


(y-y0>  =  m(x-x0) 


(2-4) 


where 

m  =  arctan  (course  line) 

Substituting  equation  (2-4)  into  equations  (2-1),  (2-2),  and  (2-3)  results  in  a  set  of 
equations  for  the  trackline  in  TD  space. 

(TDY-TDYJ  =  C^TDX-TDXq)  +  nt  (2-5) 

(TDZ-TDZq)  =  C2(TDY-TDYo)  +  n2  (2-6) 

(TDX-TDX0)  =  C3(TDZ-TDZ0)  +  n3  (2-7) 


where 


C|  =  (a2i+a22m)/(aji+ai2m) 
c2  =  (a31+a32m^a21+a12m^ 
C3  =  (an+ai2m)/(a31+a32m^ 

nl  =  ny  -  nxCj 
n2  ~  nz  ”  ny^2 
n3  =  nx  “  nz^3 


Linear  regression  can  be  used  to  fit  a  straight  line  to  TD  data  collected  along  each  of 
the  tracklines  using  a  visual  range  as  a  position  reference.  The  resultant  regression 
lines  of  TDs  are  in  the  form: 

(TDY-TDYj)  =  aj(TDX-TDX:)  (2-8) 

(TDZ-TDZj)  =  b{(TDY-TDYj)  (2-9) 

(TDX-TDXj)  =  Cj(TDZ-TDZj)  (2-10) 
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where 


TDXj,  TDYj,  TDZj  are  averages 
i  =  1,  2  (i.e.  trackline  1  and  trackline  2) 

Each  of  the  above  pairs  of  simultaneous  equations  can  be  solved  to  estimate  the 
waypoint  time  differences.  Two  estimates  for  each  of  the  waypoint  TDs  are  obtained. 
The  agreement  of  these  estimates  is  one  measure  of  quality  of  the  survey. 

The  quality  of  the  resultant  waypoint  survey  is  a  function  of  several  factors.  The 
location  of  the  ranges  is  an  important  consideration.  In  the  ideal  case,  both  sets  of 
range  markers  are  near  the  waypoint.  The  ability  of  the  operator  to  determine  when 
the  survey  vessel  is  on  the  range  decreases  as  a  function  of  distance  from  the  range 
markers.  The  confidence  bounds  cm  the  regressive  lines  are  minimum  at  the  mean 
values;  therefore  the  ideal  survey  pattern  is  an  "X"  centered  (approximately)  at  the 
waypoint.  This  pattern  is  not  always  realizable  since  there  may  be  insufficient  water 
beyond  the  waypoint,  or  range  markers  may  become  obscured  shortly  after  the 
waypoint  is  passed.  The  survey  tracklines  should  be  kept  short  to  insure  linearity  of 
the  TD  grid  over  the  survey  area.  Several  runs  are  made  on  each  trackline  to 
randomize  errors  in  positioning  the  survey  vessel  on  the  trackline. 

The  crossing  tingles  of  the  survey  lines  also  affect  the  accuracy  of  the  solution.  In 
general,  when  the  tracklines  cross  at  a  shallow  angle  cm  the  navigation  chart,  the  TD 
tracklines  will  also  cross  at  shallow  angles. 

Figures  2-1,  2-2,  and  2-3  illustrate  the  determination  of  waypoint  TDs  based  on  TDs 
measured  on  two  intersecting  ranges  on  the  St.  Marys  River.  This  is  almost  an  ideal 


FIGURE  2-1.  TDY  VS  TDX  MEASURED  ON  TWO  INTERSECTING  VISUAL 
RANGES  WITH  RESULTING  REGRESSION  LINES 
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case.  Both  sets  of  range  markers  were  near  the  waypoint  and  it  was  possible  to 
bracket  the  waypoint  in  an  "X"  survey  pattern.  The  crossing  angles  are  near  ideal 
(50-88  degrees).  The  resultant  solutions  display  excellent  agreement  as  shown  in  Table 
2-1. 

TABLE  2-1.  COMPARISON  OF  WAYPOINT  CALCULATIONS  FOR  EACH 
PAIR  OF  REGRESSION  LINES  SHOWN  IN  FIGURES  2-1,  2-2  &  2-3 

DATA  SET  TDX  TDY  TDZ 

TDX/TDY  11260.256  22332.410 

TDY/TOZ  -  22332.410  33299.094 

TDX/TDZ  11260.257  -  33299.092 


Unfortunately  not  all  channels  are  marked  by  ranges.  The  approach  in  this  ease  is  an 
interactive  one  which  utilizes  the  FEHG  algorithm  and  the  survey  officer's  judgement. 
The  first  step  is  to  survey  the  TDs  of  the  channel  features  (channel  edges, 
aids-to-navigation,  shoals,  etc.)  in  the  area  around  the  waypoint.  A  fathometer  can  be 
used  to  detect  channel  edges  and  shoals.  Buoys  and  fixed  aids  are  marked  by  circling 
or  stationing  near  them.  An  initial  estimate  is  made  of  the  waypoint  TD  from  a  simple 
prediction  program  and  the  FEHG  algorithm  can  be  used  to  plot  the  channel  features 
surveyed.  If  the  estimate  of  the  waypoint  is  incorrect,  the  location  of  the  waypoint 
with  respect  to  the  channel  features  will  appear  offset.  The  correction  which  should 
be  applied  to  correct  this  offset  in  x-y  coordinates  translates  into  a  TD  correction  for 
the  waypoint.  The  procedure  can  be  repeated  until  the  survey  officer  is  satisfied  that 
the  waypoint  is  positioned  correctly.  This  procedure  works  particularly  well  when  one 
of  the  channels  is  marked  by  a  range. 


FIGURE  2-4.  PLOT  OF  CHANNEL  FIGURE  2-5.  PLOT  OF  CHANNEL 

EDGE  FEATURES  REFERENCED  TO  EDGE  FEATURES  REFERENCED  TO 

INCORRECT  WAYPOINT  TDS  CORRECTED  WAYPOINT  TDS 
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FIGURE  2-7.  PORTION  OF 
NAVIGATION  CHART  FOR  AREA 
SHOWN  IN  FIGURE  2-6 


FIGURE  2-6.  TYPICAL  PLOT  OF 
CHANNEL  EDGE  FEATURES  DERIVED 
FROM  TD  MEASUREMENTS 


Figures  2-4  and  2-5  illustrate  the  above  approach  applied  to  a  waypoint  on  the  St. 
Marys  River  where  neither  channel  was  marked  by  a  range.  The  waypoint  in  this 
example  was  defined  as  the  intersection  of  the  centerlines  of  two  adjacent  channels. 
In  Figure  2-4  the  buoys  on  the  east  side  of  the  channel  for  the  southern  trackline  plot 
appear  to  lie  almost  on  what  should  be  the  centerline.  On  the  northern  trackline,  the 
channel  centerline  and  western  edge  coincide.  Figure  2-5  shows  the  result  of  moving 
the  waypoint  such  that  the  tracklines  are  properly  centered  in  both  of  the  channels. 

Channel  Features 

Once  the  TD  coordinates  of  a  waypoint  have  been  determined,  the  next  step  in  the 
survey  is  to  determine  the  position  offsets  of  channel  features  with  respect  to  the 
surveyed  waypoint.  A  plot  is  made  using  the  TD  data  for  the  channel  features  near  a 
waypoint  transformed  by  the  FEHG  into  planar  coordinates.  The  survey  officer  then 
digitizes  and  stores  the  relative  x-y  locations  of  channel  edges,  shoals, 
aids-to-navigation,  etc. 

Figure  2-6  illustrates  a  typical  plot  of  channel  edge  features  in  x-y  coordinates  in  the 
area  around  one  of  the  waypoints  on  the  St.  Marys  River.  Figure  2-7  is  a  portion  of  the 
navigation  chart  covering  the  same  area. 
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Daisy-chaining  Chartlets 


The  next  step  in  the  survey  process  is  to  calculate  the  position  of  each  of  the 
waypoints  based  upon  the  survey  data.  A  central,  well-established  waypoint  is  defined 
as  the  reference  point  for  the  overall  survey  area.  The  positions  of  the  waypoints 
adjacent  to  the  reference  waypoint  are  calculated  based  on  the  reference  waypoint 
position  and  TDs  and  the  surveyed  TDs  of  the  adjacent  waypoints.  Successive  waypoint 
positions  are  calculated  in  the  same  manner  from  each  adjacent  waypoint.  The 
differences  between  projected  waypoint  positions  and  the  original  positions  estimated 
from  nautical  charts  are  due  to  original  charting  errors,  chart  pickoff  errors, 
cumulative  L oran-C  grid  warp  waypoint-to-waypoint,  and  TD  survey  errors. 

Short-distance  Grid  Warps 

While  the  procedure  of  daisy-chaining  waypoints  absorbs  the  effects  of  long-distance 
grid  warps,  there  may  exist  short-distance  grid  warps  along  tracklines.  This  error  can 
be  minimized  by  introducing  trackpoints,  which  are  reference  points  located  between 
two  waypoints.  Both  the  position  and  TDs  for  a  trackpoint  can  be  estimated  from  the 
channel  feature  data. 

Survey  Equipment 

A  Time-Difference  Survey  Set  (TDSS)  was  fabricated  at  the  U.S.  Coast  Guard 
Research  and  Development  Center  to  permit  collection  and  analsis  of  Loran-C  TD 
data.  The  TDSS  consists  of  an  Austron-5000  monitor  receiver  interfaced  to  a  Hewlett 
Packard  9845  calculator.  The  TDSS  provides  real  time  data  display,  processing,  and 
storage  during  the  data  collection  phase,  and  provides  data  reduction  and  analysis 
during  post  mission  analysis.  During  data  collection,  three  sets  of  TDs  are  processed 
and  stored.  The  TDSS  provides  a  real  time  display  of  the  data  in  the  form  of  a  TD-TD 
plot  on  the  calculator's  CRT  display.  Cumulative  statistics  are  calculated  as  the  data 
is  collected  and  are  displayed  in  two  forms.  A  bar  graph  on  the  CRT  display  indicates 
the  confidence  bound  on  the  regression  lines  calculated.  The  survey  officer  may  also 
obtain  a  printout  of  the  cumulative  statistics  on  command.  At  the  end  of  a  data  run 
the  calculator  outputs  a  hard  copy  of  the  plot  and  a  printout  of  the  statistical  data.  If 
a  waypoint  is  being  surveyed  using  two  ranges,  the  survey  officer  may  calculate  the 
mean  waypoint  TDs  on  location.  Depending  on  the  results  of  the  calculation,  he  may 
choose  to  resurvey  one  or  more  of  the  tracklines  or  to  move  on  to  the  next  task.  TD 
data  collected  along  channel  edges  may  be  translated  to  x-y  coordinates  and  plotted  on 
the  CRT  on  location  for  an  on-scene  check  of  the  data  collected. 

A  fathometer  and  highly  maneuverable  shallow  draft  vessel  complete  the  equipment 
necessary  for  a  "visual  survey."  With  the  fathometer  and  suitable  vessel,  it  is  a 
relatively  simple  task  to  find  the  dredged  channel  boundaries  and  maneuver  along 
them.  The  shallow  draft  and  maneuverability  features  of  the  vessel  simplify  the  task 
of  marking  aids  to  navigation  and  collecting  data  on  ranges  outside  channel 
boundaries.  A  typical  survey  crew  is  comprised  of  a  survey  officer,  equipment 
operator/technician,  and  vessel  operator. 
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AUGMENTATION  TECHNIQUES 


The  line  where  a  harbor  area  begins  and  the  Coastal  Confluence  Zone  ends  is  a  fuzzy 
one.  Areas  such  as  the  approach  to  New  York  Harbor,  areas  of  Puget  Sound, 
Chesapeake  Bay,  Delaware  Bay,  etc.,  may  not  be  suitable  for  using  Visual  Grid  Survey. 
In  general,  such  areas  are  relatively  wide  open  and  there  are  a  minimum  of  aids 
marking  channels.  The  problem  is  more  to  survey  channel  separation  zones  rather  than 
well  marked  restricted  navigation  channels.  Inconsistencies  as  large  as  several 
hundred  feet  between  positions  estimated  from  navigation  charts  and  a  position 
reference  are  undetectable  by  the  mariner  in  such  situations.  In  these  areas  it  may  be 
more  efficient  to  use  a  survey  system  which  incorporates  a  position  reference.  Two 
such  systems  have  been  conceptualized. 

Inertial  Survey  System 

As  its  name  implies,  the  Inertial  Survey  System  (ISS)  uses  a  survey  grade  inertial 
navigation  system  (INS)  incorporating  special  operating  procedures  and  post  mission 
data  analysis  as  a  position  reference,  and  a  helicopter  as  the  survey  vessel.  Loran-C 
TD  data  is  collected  as  the  helicopter  hovers  over  the  waypoint.  The  helicopter  pilot 
receives  guidance  information  from  a  navigation  display  driven  by  the  INS.  Waypoint 
positions  input  to  the  INS  are  estimated  from  the  navigation  charts  for  the  area  being 
surveyed.  Position  corrections  to  account  for  inertial  drift  between  updates  are 
calculated  at  the  end  of  the  data  collection  period.  The  Loran-C  TDSS  used  for  visual 
survey  is  too  heavy  and  bulky  for  helicopter  deployment,  and  a  helicopter  TD 
measurement  package  will  have  to  be  developed  to  fill  this  need. 

The  ISS  using  a  helicopter  as  the  survey  platform  has  the  advantage  of  being  able  to 
cover  a  large  area  in  a  short  time.  The  biggest  disadvantage  is  cost  to  lease  and  the 
availability  of  inertial  survey  services.  Geodetic  control  points  must  also  be 
established  as  update  points  for  the  ISS.  This  adds  to  the  total  survey  costs  and  time 
to  implement. 

Microwave  Survey  Reference 

The  microwave  survey  system  conceptually  consists  of  the  Loran-C  TD  measurement 
system  used  for  the  visual  survey  system  with  an  added  position  reference  system 
input.  As  in  the  ISS  approach,  waypoint  positions  are  estimated  from  navigation 
charts.  Instead  of  hovering  over  the  waypoint,  the  survey  vessel  records  TD  and 
position  data  while  maneuvering  in  a  cloverleaf  pattern  about  the  waypoint.  The 
position  of  the  waypoint  and  the  location  of  the  vessel  are  displayed  on  the  calculator 
CRT.  Waypoint  TDs  are  calculated  by  reflecting  the  measured  TDs  to  the  waypoint 
using  a  linear  transformation. 

TDp  =  TO  +  A  (Zp-Z)  (2-11) 
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where 


TDp  is  the  estimated  waypoint  TD  vector 
TD  is  the  measured  TD  vector 
Zp  is  the  waypoint  position  vector 
Z  is  the  position  vector  at  the  measurement  point 

A  is  a  gradient  matrix  which  is  a  function  of  the  positions  of  the  Loran-C 
transmitting  stations  and  the  waypoint. 

Individual  samples  are  averaged  until  the  confidence  bound  of  the  mean  is  within  some 
preset  tolerance  level. 

Data  is  also  collected  along  traeklines  between  waypoints  and  perpendicular  to  the 
tracklines  using  the  positioning  system  as  a  reference.  This  data  is  used  to  determine 
the  presence  of  grid  warp  and  to  bound  the  navigation  errors. 

SUMMARY 

The  Visual  Survey  Method  is  the  ultimate  approach  for  Loran-C  TD  grid  survey  in 
restricted  waterways.  This  approach  ties  together  the  Loran-C  TD  grid  and  the  world 
of  the  mariner.  The  result  is  a  set  of  waypoint  position  and  TD  coordinates  and 
relative  position  coordinates  for  harbor  features  such  as  channel  edges  and 
aids-to-navigation  which  are  incorporated  into  chartlets  for  Loran-C  user  equipment. 
The  position  coordinates  for  these  channel  features  are  calculated  from  TD 
measurements  based  on  waypoint  TD  and  position  coordinates.  The  waypoint  positions 
are  calculated  in  a  daisy-chain  fashion  from  a  central  waypoint  based  on  the  difference 
in  TDs  between  adjacent  waypoints.  This  procedure  eliminates  charting  and  chart 
pick-off  error  and  absorbs  long  range  grid  warp.  The  remaining  error  sources  are  due 
to  errors  in  the  TDs  for  the  waypoints  and  channel  features  and  local  grid  warp.  The 
effects  of  local  grid  warp  are  minimized  by  establishing  trackpoint(s)  as  necessary 
between  waypoints. 

In  some  areas  it  may  not  be  feasible  to  apply  the  visual  technique.  Two  approaches 
have  been  conceptualized  to  survey  these  areas.  At  this  time  it  is  not  clear  whether 
these  approaches  will  be  required.  The  accuracy  of  relatively  simple  propagation 
models  may  be  more  than  adequate  to  calculate  waypoint  TD  coordinates  in  these 
cases. 
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Section  3 


GRID  STABILITY  ANALYSIS 


INTRODUCTION 

Precise  position  determination  with  Loran-C  depends  on  the  existence  of  a  stable, 
repeatable  time  difference  (TD)  grid.  Grid  stability  analysis,  together  with  the 
development  of  user  equipment  and  TD  surveying  techniques,  constitute  the  Coast 
Guard's  three-pronged  approach  to  harbor  and  harbor  entrance  (HHE)  navigation.  The 
Coast  Guard's  signal  stability  measurement  program  has  concentrated  on  the  St.  Marys 
River  mini-chain. 

Data  was  first  collected  over  the  mini-chain's  service  area  between  fall  1977  and 
spring  1978  as  input  to  the  design  and  calibration  of  a  grid  prediction  algorithm.6  TD 
variations  at  several  sites  were  unusually  large,  with  an  average  of  300  nanoseconds 
and  a  maximum  of  850  nanoseconds.7  These  instabilities,  if  true,  could  render  the 
mini-chain  unusable  for  precision  navigation  on  the  St.  Marys  River.  Although 
reasonable  care  was  taken  in  performing  the  measurements,  the  data  remained  suspect 
due  to  the  equipment  and  methodology  used.  A  careful  recollection  effort  was  deemed 
necessary  to  verify  the  chain's  performance  and  to  definitively  judge  its  navigation 
capabilities. 

PREPARATIONS 

Preparations  for  a  one  year  in-depth  study  of  the  mini-chain's  stability  were  completed 
in  early  May  1979.  The  major  tasks  involved  the  creation  of  three  fixed  monitor  sites 
plus  the  improvement  of  equipment  and  staffing  at  the  System  Area  Monitor  (SAM) 
station.  Figure  3-1  shows  the  layout  of  the  chain. 

Fixed  Monitor  Sites 

The  new  monitor  sites  are  on  the  south  end  of  the  river  at  DeTour  Village,  mid-river  at 
Dunbar  Forest,  and  at  Point  Iroquois  near  the  river's  northern  end.  Each  site  has  a 
35-foot  whip  antenna  with  a  groundplane  and  multicoupler.  DeTour  (Figure  3-2)  and 
Dunbar  are  equipped  with  a  Magnavox  AN/BRN-5  receiver  plus  an  Interaav  LC-204 
receiver.  Point  Iroquois  (Figure  3-3)  has  two  LC-204  receivers.  This  setup  enables 
each  site  to  monitor  all  three  baselines. 

System  Area  Monitor 

Figure  3-4  shows  the  SAM  site,  which  has  been  converted  to  a  Coast  Guard  standard 
suite.  Control  had  previously  been  achieved  using  an  Intemav  Model  303  monitor 
receiver,  a  short  whip  antenna,  and  an  early  prototype  CALOC  (Calculator  Assisted 
Loran  Controller)  system.  Since  4  May  1979,  control  has  been  performed  using  the 
35-foot  whip  antenna  and  Austron  5000-CALOC  system  that  is  becoming  standard 


FIGURE  3-1.  LAYOUT  OF  MINI-CHAIN  WITH  RECONFIGURED  M-X 


equipment  for  long-baseline  Loran-C  chains.  The  mini-chain's  group  repetition  interval 
was  slowed  from  4930  to  5930  to  satisfy  the  Austron  5000's  processing  time 
requirements,  and  has  since  been  changed  to  6980  to  preclude  interference  from  the 
Canadian  East  Coast  Chain. 

Chain  Reconfiguration 

The  final  step  in  preparation  for  the  stability  study  was  accomplished  on  9  May  1979 
when  the  Master  and  Xray  station  designations  were  interchanged.  The  switch  was 
made  to  improve  Master  station  availability.  Prime  power  failures  have  been  more 
frequent  at  the  Gordon  Lake,  Canada  site  than  at  the  Pickford,  Michigan  site.  Power 
failures  are  significant  since  the  transmitters  have  no  emergency  generators  for 
backup.  Also,  repairs  can  be  effected  more  quickly  at  the  Pickford  site  due  to  the 
shorter  travel  times  involved. 
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FIGURE  3-4.  SYSTEM  AREA  MONITOR 
PRELIMINARY  STABILITY  ANALYSIS 

The  basic  analysis  relies  upon  one  hour  averages,  termed  the  System  Samples, 
collected  at  midday  and  midnight.  The  initial  analysis  period  includes  not  only 
environmental  changes,  but  also  the  above-mentioned  refinements  in  the  monitor  and 
control  equipment.  The  data  cannot  be  considered  the  final  measure  of  the 
mini-chain's  performance,  but  are  nonetheless  informative.  Figure  3-5  shows  the 
midday  TD  fluctuations  for  one  baseline  at  the  three  monitor  sites.  The  worst  case 
variation  is  approximately  180  nanoseconds. 


Variations  Model 

The  TD  data  has  been  analyzed  by  a  rather  simple  yet  elegant  model: 
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FIGURE  3-5.  MASTER-ZULU  MIDDAY  TIME  DIFFERENCE  VARIATIONS 

where  Z  is  the  observation  matrix  containing  the  three  TD  records  at  each  of  the  three 
monitor  sites,  A  is  the  transformation  or  geometry  matrix  describing  the  position  of 
each  site  with  respect  to  the  baselines  and  the  control  station,  ATD  is  a  uniform 
change  in  propagation  velocity,  C  is  the  common  error  matrix  containing  variations 
seen  only  by  and  corrected  by  the  control  station;  and  e  is  the  matrix  containing  local 
errors,  unique  to  each  site.  The  common  variations  are  erroneous  in  the  sense  that  no 
other  site  detected  the  changes  that  were  compensated  by  the  resulting  transmitter 
timing  adjustments.  Any  non-uniform  propagation  effects  will  be  included  in  the  local 
error  terms.  The  least  squares  solution  to  the  model  iss 
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TABLE  3-1.  COMMON  CONTROL  ERRORS 


Nanoseconds 

(RMS) 

_ CX  Cy  C2 

22  Feb-8  May  1979  46  18  28 

11  May-24  Sep  1979  7  7  6 


TABLE  3-2.  TYPICAL  UNIFORM 
VELOCITY  TIME  DIFFERENCE 
CHANGES  (NANOSECONDS) 


M-X  M-Y 

M-Z 

Gros  Cap 
(North  End) 

40  -88 

125 

Six  Mile  Point 
(Near  SAM) 

-1 

-2 

9 

DeTour  Light 
(South  End) 

48 

-95 

132 

Table  3-1  lists  the  common  control  errors  extracted  from  the  data  for  the  periods 
before  and  after  the  chain  reconfiguration  and  monitor  station  upgrade.  Considerable 
improvement  can  be  seen  in  the  second  set  of  numbers. 

Figure  3-6  shows  the  uniform  velocity  of  propagation  effects  extracted  from  the  data, 
in  units  of  nanoseconds  per  kilometer  of  difference  distance  (differential 
Master-Secondary  distance  of  point  of  interest  vs  SAM).  The  propagation  component  is 
extracted  as  an  estimate,  based  upon  how  each  site  would  react  to  the  postulated 
uniform  propagation  change.  Such  a  change  is  expected  to  be  primarily  caused  by 
temperature  affecting  the  vertical  lapse  rate  of  the  index  of  refraction8*9. 

The  peak-to-peak  propagation  velocity  change  extracted  from  the  data  to  date  is 
approximately  two  nanoseconds  per  kilometer.  The  corresponding  time  difference 
variations  are  obtained  by  multiplying  the  plotted  values  by  the  difference  distances 
for  each  site.  (The  difference  distance  is  zero  for  a  site  on  the  same  hyperbolic  TD 
line  as  the  SAM  and  increases  with  distance  from  this  line.)  Typical  values  are  listed 
in  Table  3-2. 

The  uniform  propagation  changes  are  interesting  in  that  they  represent  a  fundamental 
performance  limit  in  the  system.  The  significant  point  is  that  the  information  provides 
a  bound  on  a  physical  effect  which  cannot  be  improved  upon  through  equipment  or 
procedural  changes  at  the  SAM,  but  which  might  be  improved  with  an  "altimeter 
correction"  applied  within  the  user  equipment. 


FIGURE  3-6.  UNIFORM  PROPAGATION  VELOCITY  CHANGES 
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FIGURE  3-7.  WINTER  PROPAGATION  VELOCITY  CHANGES 


Figure  3-7  shows  the  mid-February  to  mid-April  propagation  velocity  changes, 
superimposed  with  a  plot  of  temperature  recorded  at  the  SAM.  The  correlation 
coefficient  for  these  plots  is  approximately  0.7.  This  relatively  high  correlation 
confirms  the  expected  temperature  dependency.  The  correlation  diminishes  rapidly 
during  late  April  and  early  May,  after  which  there  is  negligible  correlation  with 
temperature.  This  diminishing  correlation  is  expected,  since  the  temperature 
dependency  of  the  vertical  lapse  rate  is  predominantly  a  winter  effect.8*9 

Error  Budget  Analysis 

Knowing  the  nature  and  source  of  the  grid  variations  is  interesting,  but  of  no 
immediate  value  to  the  user.  However,  having  determined  the  expected  TD 
fluctuations,  it  is  now  possible  to  calculate  the  positioning  errors  that  would  result. 
Figure  3-8  is  the  plot  of  an  error  budget  analysis  for  a  series  of  points  along  the  St. 
Marys  River.  A  user's  time  difference  error  has  been  allocated  to  three  sources: 

1.  Uniform  Propagation  Velocity  Changes.  The  maximum  variation  experienced 
to  date  (two  nanoseconds  per  kilometer  of  baseline  length)  is  used,  recognizing  that  the 
river  might  be  surveyed  during  one  extreme  and  navigated  during  the  opposite 
extreme.  The  narrower  portions  of  the  navigation  channel  are  nearer  the  SAM  and  will 
experience  the  smallest  variations,  as  can  be  seen  from  Table  3-2. 

2.  Control  Errors.  The  "erroneous"  control  station  adjustments  are  applied  as 
normally  distributed  random  numbers  with  a  zero  mean  and  a  standard  deviation  equal 
to  the  smaller  RMS  values  from  Table  3-1. 

3.  Receiver  Offsets.  Assuming  the  use  of  a  survey-quality  receiver,  IS 
nanoseconds  is  allowed  for  the  offset  between  a  user's  receiver  and  the  calibration 
receiver  which  surveyed  the  river.  The  sign  of  the  offset  is  allowed  to  vary  randomly 
for  each  of  the  three  TDs. 

Radial  error  values  we  determined  by  computing  a  minimum  variance  three-TD  "fix" 
according  to  the  "G-Matrix"  algorithm  developed  in  Appendix  B.  The  99  percent  curve 
was  determined  by  running  a  Monte  Carlo  simulation  and  plotting  the  resulting 
probability  histogram.  Also  plotted  is  the  half  width  of  the  navigation  channel.  It  can 
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be  seen  that  the  99  percent  error  is  well  within  this  half-channel  distance  for  the 
entire  river.  This  makes  an  interesting  comparison,  but  the  intent  is  not  to  imply  that 
the  half-channel  width  is  an  acceptable  error  in  a  radionavigation  system.  Allowances 
must  also  be  made  for  navigation  error  and  the  half-width  of  the  user's  vessel  (up  to  16 
meters). 

HARBOR  MONITOR  PROGRAM 

The  Office  of  Research  and  Development  is  responsible  for  evaluating  the 
effectiveness  of  the  Loran-C  system  in  providing  the  accuracy  necessary  to  serve  as  a 
reliable,  all-weather  radio  aid  within  the  HHE  environment.  The  desired  technical 
approach  is  to  extract  the  inherent  accuracy  of  the  system  through  improved 
understanding  of  the  bounding  physical  elements.  This  understanding  will  be  developed 
through  analysis,  model  development,  testing  and  verification. 

The  objective  of  the  Harbor  Monitor  Program  is  to  characterize  the  stability  of  the 
existing  grid  in  the  HHE  areas  of  the  United  States.  Quantification  of  the  year-round 
temporal  variations  is  required  as  input  to  an  error  budget  analysis  for  each  HHE  area 
of  interest. 

Several  sets  of  equipment  have  been  developed  previously  for  monitoring  the  operation 
of  both  the  West  Coast  chain  and  the  St.  Marys  River  mini-chain.  The  Coast  Guard 
Research  and  Development  Center  has  constructed  a  prototype  monitor  system  (Figure 
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3-9)  consisting  of  hardware  and  software 
to  allow  remote  dial-up  access  to  the 
daily  System  Samples.  These  prototypes 
are  presently  deployed  at  the  St.  Marys 
River  monitor  sites. 

Improved  versions  of  the  remote  moni¬ 
toring  hardware  will  be  installed  in  each 
harbor  area  studied.  Up  to  ten  sets  will 
be  installed  during  1980.  Two  monitoring 
sites  should  be  sufficient  to  characterize 
the  grid  stability  in  most  regions.  Boston 
Harbor,  the  Delaware  Bay  and  the 
Chesapeake  Bay  have  been  chosen  for 
the  first  installations. 

CONCLUSIONS  FIGURE  3-9.  HARBOR  MONITOR  SYSTEM 

The  preliminary  results  have  established  that  the  grid  is  reasonably  stable,  and  that 
earlier  reports  showing  much  larger  variations  were  apparently  incorrect. 

An  analysis  of  the  data  has  revealed  a  limiting  physical  effect  in  the  form  of  changes 
in  the  velocity  of  propagation,  uniform  throughout  the  chain.  These  propagation 
changes  have  been  as  large  as  130  nanoseconds  at  the  ends  of  the  river,  but  are  much 
smaller  near  the  control  station. 

An  error  budget  analysis  has  projected  three-TD  fixes  that  are  well  within  the 
half-width  of  the  river's  navigation  channels.  However,  the  addition  of  navigation 
error  and  the  half-width  of  a  typical  vessel  preclude  making  the  implication  that  the 
mini-chain  is  capable  of  providing  safe  electronic-only  navigation  in  every  section  of 
the  St.  Marys  River  at  all  times. 

PROGNOSIS  AND  FUTURE  EFFORTS 

The  mini-chain  stability  study  has  recently  been  completed  (May  1980),  and  analysis  of 
the  additional  data  collected  since  September  1979  should  enable  a  definitive 
judgement  to  be  made  on  the  navigation  capabilities  of  the  St.  Marys  River  mini-chain. 

The  geometry  of  the  Great  Lakes  Loran-C  chain  is  excellent  in  the  St.  Marys  River 
area.  Consideration  is  being  given  to  providing  precision  coverage  of  the  river  by 
supplementing  the  long-baseline  chain  with  a  single  mini-Loran-C  station  —  a  much 
more  cost  effective  approach  —  in  lieu  of  providing  an  entire  mini-chain.  Present 
plans  call  for  terminating  the  mini-chain  concept  after  the  St.  Marys  River  is  closed  to 
navigation  for  the  1980-1981  winter  season.  The  Gordon  Lake  station  will  then  be 
brought  on-air  as  a  low-powered  secondary  on  the  Great  Lakes  rate  (8970).  A  new 
experiment  will  be  developed  to  verify  the  validity  of  this  augmented  coverage 
concept.  Data  collection  for  this  new  experiment  was  started  on  the  Great  Lakes 
chain  at  several  of  the  St.  Marys  River  monitor  sites  during  the  1979-1980  winter 
season,  and  was  recently  expanded  to  include  all  sites. 
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TO:  J.  G.  Wall 

FROM:  L.  J.  Levy 


SUBJECT:  Definition  of  G  Matrix  Analysis  for  USCG  LORAN-C 

HKE  Navigation  Program 

REFERENCES:  1)  APL/JHU  Technical  Program  Plan  for  the  USCG 

LORAN-C  Harbor-Estuary  Navigation  Program, 
Enclosure  (1)  of  AD-7599 

2)  LORAN-C  HHE  Navigator,  Handwritten  Notes  by 
Cmdr.  Don  Feldman,  USCG 


The  program  plan  for  the  USCG  LORAN-C  Harbor-Harbor - 
Entrance  (HHE)  Navigator  program  in  Ref.  1  states  that  APL  is  to 
"Perform  the  'G*  matrix  analysis...".  Based  on  discussions  with 
Lt.Cdr.  Dan  Garrison,  USCG,  discussions  with  the  USCG  at  the 
December  4  meeting,  and  the  notes  in  Ref.  2,  the  G  matrix  analysis 
to  be  performed  by  APL  can  now  be  defined. 

The  objective  of  the  G  matrix  study  is  to  j  commend  to 
the  USCG  an  algorithm  to  calculate  the  G  matrix  under  the  follow¬ 
ing  three  conditions: 

a)  2TD  Case 

Given  the  location  of  the  master,  two  slave  stations, 
and  a  surveyed  point,  calculate  the  unique  G  matrix 
for  that  surveyed  point. 

b)  3TD,  Weighting  Vector  Casa 

Given  the  location  of  the  master,  3  slave  stations, 
and  a  surveyed  point  and  a  set  of  TD  weighting  para¬ 
meters,  reflecting  the  confidence  in  each  TD,  calcu¬ 
late  a  G  matrix  for  that  surveyed  point  using  a 
weighted-least-squares  criterion  on  the  TD  residual. 
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c)  3TD  Statistical  Minimization  Case 

Given  the  location  of  the  master,  3  slave  stations, 
and  a  surveyed  point,  calculate  the  G  matrix  for 
that  surveyed  point  such  that  the  statistical  mean 
square  error  in  estimating  X  and  Y  is  minimized. 

It  is  to  be  assumed  that  the  noise  level  into  the 
LORAN  receivers  is  equal  and  Gaussian.  APL  will 
provide  an  algorithm  for  computing  the  required  TD 
error  covariance  matrix  to  be  primarily  based  upon 
calculated  signal  strength  and  assumed  background 
noise  considerations. 

Additionally,  it  is  assumed  that: 

1)  All  G  matrices  will  transform  TDs  in  units  of  nano¬ 
seconds  to  position  in  the  X-Y  coordinate  system 

in  units  of  yards. 

2)  Since  it  was  agreed  that  the  generation,  tape  stor¬ 
age,  an  use  of  new  G  matrix  coefficients  could  be 
easily  accomplished  for  each  waypoint,  no  linearity 
analysis  or  sensitivity  analysis  need  be  performed 
over  tracklines  not  including  the  waypoint  used  in 
the  calculation  of  a  particular  G  matrix.  Only  a 
linearity  analysis  over  the  adjacent  tracklines  will 
be  performed  to  determine  if  an  additional  waypoint, 
with  associated  new  G  coefficient,  is  needed. 

3)  APL  will  also  provide  to  the  USCG  an  algorithm  to 
evaluate  linearity  and  an  algorithm  to  evaluate  the 
covariance  performance. 

4)  A  limited  amount  of  covariance  performance  analysis 
will  be  done  for  each  case  to  show  the  utility  of  the 
covariance  analysis  algorithm. 

5)  A  Geometrical  Dilution  of  Precision  (GDOP)  number 
will  be  calculated  for  each  waypoint. 

6)  The  algorithms  presented  to  the  USCG  will  be  in 
standard  mathematical  equation  form.  Non  optimized 
coding  of  these  algorithms  may  be  supplied  if  re¬ 
quested  by  the  USCG. 

7)  Grid  warpage  will  be  controlled  by  insertion  of  ad¬ 
ditional  trackpoints  to  be  determined  by  the  USCG. 
Decisions  as  to  the  choice  of  G  matrix  algorithm  is 
the  responsibility  of  the  USCG. 

The  APL  report,  recommending  the  algorithms  for  the  three 
above  cases,  plus  providing  the  associated  linearity  and  covariance 
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analysis  algorithms  will  be  submitted  to  the  USCG  by  January  15, 
1979.  A  follow-on  report  providing  the  linearity  analysis  of  item 
(2),  the  limited  covariance  performance  analysis  of  item  (4),  and 
the  GDOP  analysis  of  item  (5)  will  be  provided  to  the  USCG  by 
February  30,  1979. 
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G  Matrix  Analysis  for  USCG  Loran-C  HHE  Navigation 
Program 

1)  Loran-C  HHE  Navigator,  Handwritten  Notes  by 
Cmdr  Don  Feldman,  USCG 

2)  Definition  of  G  Matrix  Analysis  for  USCG  Loran-C 
HHE  Navigation  Program,  S2R-78-307,  by  L.  J.  Levy 
December  6,  1978 

3)  C-LAD  -  A  Low  Cost  Loran-C  Assist  Device;  Vol  II 
Software,  L.  M.  Marshall,  C.  R.  Edwards,  APL/JHU 
CP034B,  June  1974 


INTRODUCTION 

The  G  Matrix  technique,  as  defined  in  Reference  1,  is 
essentially  a  method  for  converting  2  or  more  measured  TDs  into 
horizontal  navigator  position.  It  yields  a  linear  estimate  refer¬ 
enced  to  the  surveyed  values  of  a  "near"  waypoint  or  track  point 
as 


2  =  Z  +  G  TtD  -  TD  "1  (1) 

-P  1_  Q  “Pj 


where 

Z_  -  surveyed  position  vector  of  point  P,  a 
P  nearby  TP  (track  point)  or  WP  (waypoint) 

^  '• 

Z^  ~  estimate  of  receiver  position  vector  at 

Q  point  Q 

TD  ~  vector  of  surveyed  time  differences  at 
P  point  P 

TD  -  vector  of  measured  time  differences  from 
Q  receiver  at  point  Q 


TO: 

FROM: 
SUBJECT : 

REFERENCES : 
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G  -  G  matrix  which  transforms  the  vector  TD 

residual  into  vector  position  residual  in 
some  optimal  fashion 

In  general,  the  G  matrix  is  a  function  of  the  surveyed 
position  of  the  WP  or  TP,  surveyed  position  of  the  master  and 
slave  stations,  an  assumed  estimation  error  minimization  criterion, 
and  possibly  of  the  expected  relative  error  characteristics  or 
quality  factors  of  TD^.  Its  validity  is  based  upon  the  assumption 
that  the  residual  TD  vector  (measured  TD  vector  minus  the  surveyed 
TD  vector)  is  linearly  related  to  the  position  residual  (position 
of  navigator  minus  surveyed  position  of  the  WP  or  TP) . 

The  significance  of  the  G  matrix  technique  is  that  the 
G  matrix  can  be  precomputed  for  each  WP  or  TP,  thereby  significantly 
reducing  the  amount  of  real  time  computation  of  the  navigator  in 
transforming  measured  TDs  into  navigator  position. 

The  objective  of  the  G  matrix  analysis  is  to  recommend 
to  the  USCG  the  algorithms  to  calculate  the  G  matrix  under  the 
three  conditions  specified  in  Reference  2.  This  report  will  present 
these  recommended  algorithms  along  with  the  associated  covariance 
and  linearity  analysis  algorithms. 

1.  General  G  Matrix  Solution 

This  section  will  define  the  navigation  problem  and  pre¬ 
sent  the  general  G  matrix  solution.  Given  an  m-vector  of  measured 
TDs,  TD  ,  at  navigator  position  Q;  the  corresponding  vector  of 
surveyed  TDs,  TD^ ,  for  a  nearby  point  P;  and  the  surveyed  position 
2-vector,  Z^,  of  the  nearby  point  P;  determine  the  position  of  the 
navigator.  From  Appendices  B  and  C,  it  is  shown  that  when  point  Q 
is  sufficiently  near  point  P, 

u  =  Ax  +  e  (2) 
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where  A  is  a  function  of  the  positions  of  the  master,  slaves,  and 
the  surveyed  point  P, 

U  =  TD  -  TD  (3) 

Q  P 


X  =  Z_  -  Z  (4) 

Q  P 


and  e  is  the  m-vector  of  TD  measurement  errors  at  point  Q.  Since 

A 

Z_  is  known,  then  an  estimate  of  x,  x,  will  yield  an  estimate  of 


by 


As  shown  in  Appendix  A,  the  general  form  of  the  estimate  of  x  is 
given  by 


x  =  G  u  (6) 


where 

G  =  <ATW  A)_1ATW  (7) 

For  the  2  TD  case,  m=2,  the  G  matrix  reduces  to 

G  =  A*1.  (8) 

Otherwise,  the  specific  value  of  the  G  matrix  is  dependent  upon 
the  weighting  matrix,  W.  If  W«I ,  the  measurements  are  equally 
weighted  and  Eqs.  {€)  and  (7)  yield  the  standard  least  squares 
estimate  of  x.  For  a  general  W,  we  have  a  weighted  least  squares 
estimate.  For  W=R-*,  where  R  is  the  covariance  of  e,  we  have  the 
minimum  variance  estimate. 
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Combining  Eqs.  (3),  (5),  and  (6)  we  obtain  the  desired 
result  of  estimating  the  navigator  position  at  point  Q  via  Eq.  (1) . 

2.  Specific  Algorithm  for  the  2  TD  Case 

Given  the  location  of  the  master  station,  (x  ,y  )  ;  the 

m  m 

slave  stations,  (x  .,y  .),  j-1,2,3;  and  the  surveyed  position  of 

s  D  SD 

the  WP  or  TP,  (xp,yp) ;  calculate  the  unique  G  matrix  for  each  2  TD 
case.  Here  (x,y)  are  the  components  in  the  planar  coordinate  frame.* 
(For  the  St.  Marys  River  system  the  coordinate  center  (i.e.,  point 
of  tangency)  is  at  the  Loran  monitor  station  (46°  28.1244*  N, 

84°  17.9343'  W) ) .  At  the  coordinate  center  x  is  east  and  y  is 
north.  For  the  2  TD  case,  G  is  given  by 


G(i, j)  = 


gxi 

gxj 

gyi 

gyj 

th 

station 

(9) 


coordinate  (a=x,y) .  The  G(i,j)  matrix  transforms  the  i  and  j  sta¬ 
tion  TDs  into  the  x  and  y  position  estimates.  Now  from  Appendix  B, 


ri  -y<vv’  *  (vyi»)2  <io) 

+  vysj>2 !  j  ■ 1,2,3  tii) 

where  r^  is  the  slant  range  from  the  surveyed  point  P  to  the 
master  station,  r2j  is  the  slant  range  from  the  surveyed  point  P 
to  the  jth  slave  station.  The  gradient  matrix  values  are 


*  Planar  coordinate  frame  defined  in  Reference  3,  pp.  30-35,  141-144; 
x*p  sin¥,  y«p  cos¥,  where  p*=geodetic  arc  length  from  origin  to 
desired  point,  Sheading  angle  of  geodetic  arc. 
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i  j 

-  1/2,3 

(12) 

t  j 

=  1,2,3 

(13) 

where 


Then 


a^x  is  the  derivative  of  the  j^“  TD  with  respect  to  x^, 
a_.y  is  the  derivative  of  the  jth  TD  with  respect  to  yp. 


v  is  the  propagation  velocity. 


A (i , j)  -a.  a.  -a.  a. 

ix  }y  iy  jx 


for  i=l  j=2 ,  i=l  j  =3 ,  i=2  j  =  3 


and  for  G(i,j) 


ajy/Mi,  j) 


-a.y/A(i,j) 


-a^/ACi,  j) 


for  i=»l  j=2 
i-1  j=3 
i*2  j-3 


9yi  = 


Thus,  3  separate  and  distinct  G  matrices  are  possible  for  the  2  TD 
case.  Note,  for  instance,  that  gxl  in  G(l,2)  does  not  equal  gxl 
in  G(l,3).  Finally,  note  that  if  x  and  y  are  in  units  of  yards, 
v  in  units  of  yards/nanoseconds,  then  the  elements  of  the  G  matrix 
are  in  units  of  yards/nanoseconds.  Equations  (9)  through  (15) 
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comprise  the  recommended  algorithm  for  calculating  the  three  G 
matrices  for  the  2  TD  case. 


3.  Specific  Algorithm  for  the  3  TD,  Weighting  Vector  Case 

Given  the  location  of  the  master,  3  slave  stations,  a 
surveyed  position  at  point  P,  and  a  set  of  TD  weighting  parameters. 


(16) 


reflecting  the  confidence  (i.e.,  large  implies  high  confidence 
in  the  ith  TD)  in  each  TD,  calculate  the  G  matrix  for  the  surveyed 
point  P  using  a  weighted- leas t-squares  criterion. 

First  of  all,  it  must  be  noted  that  the  weighting  para¬ 
meters  do  not  need  to  be  normalized  (i.e.,  w^+w2+w3=l)  although  it 
may  be  meaningful  from  an  interpretation  standpoint.  This  can  be 
seen  by  substituting  aW  for  W  in  Eq.  (7)  where  a  is  a  scalar.  The 
value  of  G  remains  the  same. 

For  this  case  we  set  the  weighting  matrix  by 


w 


1 


W 


0 


0 


(17) 


We  see  by  Appendix  A  that  this  will  result  in  weighting  the  square 
of  each  TD  residual  according  to  our  respective  confidence  in  each 
TD  in  the  weighted- squared- residuals  function  that  is  minimized. 
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The  resulting  least  squares  estimate  is  given  by  the  G  matrix  in 
Eq.  (7)  using  Eq.  (17).  The  specific  result  using  Eq.  (17)  in 
Appendix  C  and  taking  into  account  the  zeroes  in  W, 


b*j 


V 


a .  w . 
jx  3 


a .  w . 
DY  3 


j=l,2,3 
j-l»  2,3 


:  =  Z  b  .a. 

XX  j=1  X]  ]X 


;  =  Z  b  .  a . 

xy  j=1  X3  ]y 


2  =  Z  b  .  a.  «*  c 

yx  ^ml  Y3  jx  xy 


3yy  =  byj  ajy 


■  c  c  -  c 
yy  xx  xy 


-1 


XX 


-1 


xy 


c  /A 

yy 


-c  /A 
xy 


-1 


-1 


yx 


-l 


yy 


xy 


c  /A 

XX 


(18) 

(19) 


(20) 


(21) 
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then 


g  .  =  c  "1b.+c  -1  b  .  j=l,2,3  ^ 

’xj  xx  xj  xy  yj 


g  .  *  c  ^b.+c  ^  b  .  j=l,2,3 

yyj  yx  xj  yy  yj  J 


(22) 


Note  that  the  a.  ,  a.  ,  g  . ,  and  g  .  are  defined  in  Section  2.  The 

3*'  3y'  *x;j'  yyj 

resulting  G  matrix  for  this  case  is 


G(l, 2 , 3)  = 


gxl  g 


L  y1  yy2 


x2 

gx3 

y2 

gy3 

(23) 


Equations  (10)-(13),  (18) — (23)  comprise  the  recommended  algorithm 
for  calculating  the  G  matrix  for  the  3  TD  weighting  vector  case. 


4.  Specific  Algorithm  for  the  3  TD  Statistical  Minimization  Case 

Given  the  location  of  the  master,  3  slave  stations,  and  a 
surveyed  point  P,  calculate  the  G  matrix  for  that  point  such  that 
the  variance  of  the  error  in  estimating  the  navigator  position  is 
minimized.  It  is  to  be  assumed  that  the  noise  into  the  4  tracking 
loops,  measuring  receipt  times  of  the  3  slave  and  master  pulses,  are 
equal  and  Gaussian. 

From  Appendix  A,  the  minimum  variance  criterion  results  in 
setting  W*R_1  in  Eq.  (7) .  An  accurate  determination  of  R  requires 
specific  knowledge  of  the  signal-to-noise  ratios  at  the  receiver 
antennae  (which  are  receiver  configuration  independent) ,  and  many 
additional  parameters  that  are  specific  to  the  receiver  configura¬ 
tion.  Thus,  each  different  receiver  configuration  would  require  a 
detailed  error  analysis  to  ultimately  yield  a  specific  G  matrix. 

We  believe  this  to  be  impractical  and  will  strive  to  obtain  a  G 
matrix  for  this  case  that  is  configuration  independent. 
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The  constraint  that  the  noise  into  the  4  tracking  loops 
are  equal  can  be  justified  under  the  following  assumptions: 

a)  The  dominant  noise  into  the  tracking  loop  is 
external  noise  (receiver  noise  is  negligible) . 

b)  RF  receiver  input  bandwidths  are  equal  (if  the 
RF  gain  from  antenna  input  to  the  tracking  loop 
differs  between  receivers,  the  signal-to-noise 
ratios  for  the  external  noise  remain  the  same, 
yielding  the  same  relative  noise  performance  be¬ 
tween  receivers) . 


Also,  it  is  assumed  that  the  relative  noise  performance 
between  tracking  filters  from  the  tracker  filter  input  to  the  input 
to  the  G  matrix  transformation  remains  invariant.  This  can  be 
satisfied  if: 

1)  all  tracking  filters  are  equal, 

2)  quantization  noise  in  the  counting  process  is 
negligible, 

3)  any  further  smoothing  beyond  the  counting  process 
is  equivalent  between  receivers. 


Under  assumptions  (a)-(c),  (l)-(3),  the  G  matrix  is  receiver  con¬ 
figuration  invariant  and  dependent  only  on  the  effective  transmitted 
power  of  each  slave  station,  Pi#  i=l,2,3,  the  effective  transmitted 
power  of  the  master  station,  and  the  slant  ranges,  r2j;  j®1/2*3  and 
r^.  From  Appendix  D,  the  calculation  is 


2  2 
r2i  rl 

*ii  -  <-r-  +  '  i*1'2'3 

11  1  *i  m 


Rij  =  kl  (’7_)  ifil 

xa 


(24) 


where  k^  is  arbitrarily  set  so  that  the  above  elements  of  R  are 
roughly  around  unity.  (It  is  shown  in  Appendix  D  that  only  the 
relative  values  of  R^j  will  affect  the  G  matrix.) 


B-9 


THE  JOHNS  HOPKINS  ONIVfcRSi  I V 
APPLIED  PHYSICS  LABORATORY 

i  Ail  HI  (  M AH/lAHb 


S2R-79-008 
Page  -10- 


Then 


with 


-1 


„  -1  _  —  1  „  -1 

R11  R12  R13 


R  -1  R  -1  R  -1 

R12  R22  K23 


R  "I  R  "I  R  “I 

k13  r23  33 


(25) 


*  R11R22R33  +  2  R12R23R13 
"  R132r22  "  R122r33  “  R11R23 


R11  “  ^R22R33  ~  R23 

R12  1  *  (R13R23  "  R33R12)/A  =  R21 


r13-! 


-1 


22 


^R12R23  "  R22R13)/,A  ~  R 

(R11R33  “  R132)/A 


-1 


31 


-1 


23 


-1 


'33 


(R12Ri3  -  RnR23)/A  R 


(R11R22  “  R12 


-1 


32 


(26) 


From  Appendix  C,  substituting  M=R  \ 

bxj  "  aix  Rij  ?  j  *  1,2,3 

byj  ’  aiy  *ifl  J  "  X'2'3 


(27) 
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The  rest  of  the  required  equations  are  given  by  Eqs.  (20) -(23). 

Equations  (10)-(13)f  (24),  (26),  (27),  ( 20)  —  ( 23)  comprise 
the  recommended  algorithm  for  calculating  the  G  matrix  for  the  case 
considered  in  this  section. 

5.  Linearity  and  Covariance  Analysis  Algorithms 

From  Appendix  A,  the  general  form  of  the  error  covariance 
T 

is  C=GRG  ,  where  G  is  any  G  matrix  of  the  form  in  Eq.  (7)  with  any 
W,  and  R  is  the  actual  covariance  of  TD  measurement  error  vector 
(not  necessarily  the  R  matrix  computed  in  Section  4) .  Define 


Ax  = 

A 

XQ  “ 

XQ 

li 

A 

yq 

then 

— 

<Ax^ 

> 

<Ax  Ay> 

c  = 

~  7 

<  A~x 

Ay> 

<Ay  > 

For 

Xi 

Rtj‘ 

R(i,j)  = 

-Rij 

Ui- 

UJ. 

1 _ 

with  G(i,j) 

given  by  Eq. 

(9) 

,  then 

C(i,j)  =  G(i,j)R(i,j)G(i,j)T  (28) 

for  any  pair  of  TDs,  i,j. 
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For  any  3  TD  case,  set 


with  G  given  by  Eq.  (23) ,  then  C  =  GRGT. 

The  linearity  of  G,  generated  for  surveyed  point  P, 
evaluated  at  point  Q  as  follows: 

1)  Use  Gp  computed  for  particular  2  TD  case,  along 
with  surveyed  TDs  at  point  Q  and  point  P,  and 
surveyed  position,  Zp,  in  Eq.  (1)  to  get  Z^. 

2)  Linearity  error  at  point  Q  relative  to  point  P 
is  given  by  differencing  the  surveyed  position 

with  Zq. 

3)  Repeat  steps  (1)  and  (2)  for  each  2  TD  case  and 
single  3  TD  case. 
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Appendix  A 
LINEAR  ESTIMATION 

(Derivations  of  Formulae  Due  to  S.  T.  Haywood) 


A  set  of  m  linear  equations  in  n  <  n  unknowns  will  be 
represented  by 

(A.l)  u  =  Ax  +  e 

in  which  we  define: 

A  mxn  matrix  of  rank  n 
e  nxl  noise  vector 

x  mxl  vector  of  true  values  of  the  state  variables 
u  nxl  vector  of  observations  or  measurements 

In  (A.l)  u  and  A  are  presumed  known;  e  and  x  are  unknowns. 

The  problem  to  be  considered  next  is  that  of  obtaining 
an  estimate  of  x.  Let 

(A.  2)  d  *  Ax  -  u 

where 


d  nxl  vector  of  residuals 
x  mxl  vector  representing  an  estimate  of  x 

In  classical  least  squares  estimation  an  estimate  of  x  is  obtained 
which  minimizes  the  scalar 
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dTd  =  jd|  =6  +6  +  ...  +  6  . 

12  n 

We  shall  now  derive  the  weighted  least  squares  estimate 

—  T  T 

of  x  which  minimizes  d  Wd,  where  W=W  >  0  is  the  weighting  matrix. 

Ordinary  least  squares  estimate  will  be  obtained  as  the  special 

case,  W=I . 

From  (A. 2)  we  obtain 
dTWd  =  (Ax-u)TW(Ax-u)  =  ( xTAT-uT) W (Ax-u) 

W  Ip  fTl  ip  ip  Ip 

=  x  A  WAx  -  x  A  Wu  -  u  WAX  +  uWu 

Ip  ip  ip  ip  Ip  Ip  ip 

=  x  (A  WA)x  +  u  Wu  -  x  A  Wu  -  u  W  Ax 
=  xt(atwa)x  +  uTwu  -  2xtatwu 

Let  s=s(x)»dTWd.  Then  s(x+h)  =  (d+Ah)TW(d+Ah)  and 

s (x+h) -s (x)  =  dTWd+hTATWd+dTWAh+hTAWAh-dTWd 
=  hTATWAh+2hT(ATWAx-ATWu) 

The  first  term  on  the  right  above  is  always  non-negative.  The 
second  term  vanishes  if 

(A. 3)  x  -  (ATWA)"1ATWu 


Thus  if  x  has  the  value  given  by  (A. 3)  s(x+h)-s(x)  0  for  any  mxl 

vector  h.  Therefore  x  is  the  desired  weighted  least  squares  esti¬ 
mate  of  x.  The  ordinary  (unweighted)  least  squares  solution  is 
given  by 

(A.  4)  x  «  (ATA)'1ATu. 
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We  next  consider  the  statistical  properties  of  the 
noise  vector  e  which  enters  into  the  determination  of  the  co- 
variance  matrix  associated  with  the  error  in  the  estimate  of 
x.  We  shall  assume  <e>=0  and  <ee  >=R>0.  In  other  words,  we 
assume  the  noise  in  the  measurements  has  zero  mean  and  covari¬ 
ance  matrix  denoted  by  R. 

Now,  from  (A. 3), 

X  =  (ATWA)_1ATWu  =  (ATWA)“1ATWAx-(ATWA)“1ATWe 

=  x-(ATWA)_1ATWe. 


Therefore 

(A. 5)  <x>  =  x. 

Equation  (A. 5)  asserts  that  x  is  an  unbiased  estimate  of  x. 

The  covariance  matrix  of  the  error  in  the  estimate  x 
is  denoted  by  C  and  given  by 

(A. 6)  C  *  <(x— x) (x— x)T>  =  (ATWA) ”1ATWRWA(ATWA) -1. 


We  next  consider  the  problem  of  finding  the  estimate  x 
which  makes  C  a  minimum. 

Let  x  *  Bu  *  BAx  +  Be 

Then  <x>  «  B<u>  ■  BAx  *  x 

Therefore  BA  «  I 

x  -  x  ■  Bu-x  *  BAx  +  Be-x  ■  Be 
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C  =  C (B)  =  < (x-x)  (x-x) T>  =  BRBT 


Next,  let  H  be  any  matrix  having  the  same  dimensions  as  B, 


I  =  (B+H) A  =  I  +  HA 


Therefore  HA  =  0 


Then  C (B+H) -C (B)  = 


(B+H) R(BT+HT)  -  BRBT 


brbt  +  brht  +hrbt  +  HRHT  -  BRBT 


T  -1 

Let  B  =  MAR 


Then  BRHT  +  HRBT 


T  T  T 

HRH  +  BRH  +  HRB 


T  -1  T  -IT 

MAR  XRH  +  HRR  AM 


mm  T 

MA  H  +  HAM 


Since  HA  =  0,  C(B+H)-C(B)  =  HRH  >  0 


I  =  BA  =  MATR*'1A 


Therefore  M  =  (ATR-1A)-1 
Therefore  B  ■  (ATR-1A)  "^R-1 


and,  the  estimate,  x,  which  minimizes  C  is  given  by 

(A.  7)  x  =  (ATr“1A)"1ATr"1u. 

The  minimized  covariance  matrix  is  given  by 


(A.  8)  Cmin  =  (ATr"1A)"1 


Then 
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The  estimate  given  by  (A. 7)  is  known  as  the  minimum 
variance  estimate.  It  may  be  noted  that  the  minimum  variance 
estimate  is  obtained  as  the  weighted  least  squares  estimate 
when  W  *  R~*. 

The  general  expression  for  the  covariance  matrix  C 
is  given  by 

(A. 9)  C  =  BRBT 

in  which  B  can  be  any  nxm  matrix  of  the  form 

T  -IT 
B  =  (A  WA)  A  W 

In  the  special  case  where  m®n  (number  of  equations  equal 
to  number  of  unknowns)  and  the  matrix  A  in  (A. 3)  and  (A. 4)  is  in¬ 
vertible,  then 

(ATWA)"1ATW  =  A"1(ATW)~1(ATW)  =  A-1 
Accordingly  (A. 3)  and  (A. 4)  reduce  to 
(A. 10)  x  =  A~*u 

and  (A. 9)  reduces  to 
(A.  11)  C  =  A_1R(a"’1)T 


when  m=n. 
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Appendix  B 
2  TD  Case 

In  Fig.  B-l  the  master  station  is  at  M,  a  slave  station 
at  S j  and  an  arbitrary  way  point  at  P. 


M 


P 


Fig.  B-l 

Since  we  are  considering  2  slave  stations  the  index  j  assumes 
the  values  1  and  2. 

The  coordinates  of  M  will  be  denoted  by  (xmrym) «  p  ky 

(x  ,y  )  ;  S.  by  (x  . ,y  .) .  The  three  slant  ranges  shown  in  Fig. 

P  P  3  s3  SJ 

B-l  are  given  by 
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J(x  -xj2  + 
y  P  m 

<  Vym>  2 

* 

N) 

u. 

II 

I  (x  -x  . )  2  + 
1  P  S] 

W2 

,  j  =  1,2 

11 

•n 

m 

H 

_ 1 

1 (x  .-X  ) 2  + 

I  8:1  m 

(ys j-ym>  2 

,  j  =  1,2 

The  time  differences  in  the  signal  receipt  times  at  P 
are  denoted  by  and  t2  an<*  given  by 


(B.  2) 


T1  =  (r31+r2rrl)/v  +  6t1 
t2  =  (r32+r22‘rl)/v  +  +6t2 


where  the  Sts  denote  the  absolute  emission  delays  of  the  slave 
stations  and  v  denotes  the  speed  of  signal  propagation. 

The  derivatives  of  and  T2  with  respect  to  xp  and  y 
appear  as  elements  in  a  2x2  A-matrix.  These  are  given  in 


(B.3) 


^1  1  K-xsi 

*11  dxp  v  I  r21  rx  j 


12 


^ 1  _  1 
dyp  ~  v 


dx, 

l21  =  dT 


a 


22 


dT2 


1 

v 


1 

v 


^p~^sl  ^p  ^m 
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Let  Ax  and  Ay  denote  the  coordinates  of  an  arbitrary 
point  Q  (i.e.,  navigator)  relative  to  the  way  point  P  as  shown 
in  Fig.  B-2. 


Q 


Ay 


Then  the  time  differences  observed  at  Q  can  be  expressed  by 


v 

"Tl“ 

all  a12 

“Ax  ” 

"el" 

(B.  4) 

= 

+ 

+ 

_T2_ 

_T2_ 

_ a21  a22  _ 

_Ay  _ 

_e2_ 

in  which  e1  and  e2  represent  noise  terms  present  in  the  observa¬ 
tions  of  x i  and  x2  made  at  Q. 

In  (B. 4)  let 


- 1 

_ 1 

i 

H 

_ 1 

T1Q  "  T1P 

- 

s 

LTd 

% 

-T2- 

1 

H 

M 

o 

1 

H 

to 

*d 

Q  P 
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A  » 


all 

_a21 


A 


X  ■ 


Ax 

Ay 


Then  (B.4)  is  seen  to  have  the  linear  form 

u  «  Ax  +  e 

of  Appendix  A  in  which  u  and  A  are  known  and  x  and  e  are  unknowns. 
In  Appendix  A  it  is  shown  that 

(B.5)  x  «  A-1u 

is  the  least  squares  estimate  of  x  when  m=n=2. 

We  shall  next  proceed  to  write  explicit  formulas  for  Ax 
and  Ay.  The  G  matrix  is  defined  by 
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(B.6)  X  -  Gu 

Therefore, 

(B. 7)  G  =  A-1 

where 

A  =  an&22  "  a12a21 

an  '  a22/4  *  ’11 
a12  *  *  a12/4  *  912 
a2l  =  “  a21/A  =  921 
a22  “  all/A  *  g22 

and  where 

“gU  *12 

G  = 

_g21  922 

Therefore, 

(B.8)  Ax  «  *  g12^T2Q'T2P) 

Ay  -  921(t1q“T1P)  +  g22*T2Q“T2P* 
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Appendix  C 
3  TD  Case 

(based  on  definitions  and  assumption  stated  in  Appendix  B) 


In  the  case  where  3  TD  signals  are  employed  the  slant 
range  equations  are  given  in 


(C.l) 


r2j 


r3j 


fv1'"’2  + 

fVx=  i 


(Yp-Tm> 


)2  ♦  W2  . 


=  fX sj- V  2  +  (ySj-ym>  2  «  S'1'2'3 


The  three  time  delays  are 


(C.  2) 


(C.  3) 


Tj  “  (r3j+r2j_rl)/^V  +  6Tj  '  3=lf2,3 


The  six  elements  of  the  3x2  A-matrix  are 


1 

a  .  „  =  — 

x  -x  .  x-x 

p  SI  .1  rn 

jl  V 

L  r2: 

"l  J 

a ,,  =  i 

Vysj  yp"yml 

3  j  2  v 

— 

T-» 

CM 

Lc 

_ 1 

r>  J 

,  j=if 2 , 3 


j-1,2,3 


The  linear  expansion  of  the  three  measured  time  dif¬ 
ferences  made  at  an  arbitrary  point  Q  relative  to  a  way  station 
at  P  is  given  by 
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Equation  (C.4)  will  be  recognized  as  the  linear  form 
u  =  Ax  +  e 


of  Appendix  A  in  which 
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In  Appendix  A  it  is  shown  that  the  weighted  least 
squares  estimate  of  x  and  the  minimum  variance  estimate  have 
the  same  algebraic  form  which  we  shall  represent  here  by 

(C.5)  x  =  (ATMA)_1ATMu 

To  obtain  the  weighted  least  squares  estimate  we  set  M=W,  the 

least  squares  weighting  matrix.  The  minimum  variance  estimate 

—  -It 

of  x  is  obtained  when  M=R  where  R  =  <ee  >. 

We  shall  next  obtain  explicit  formulas  for  the  solu¬ 
tion  of  (C.4)  in  the  form  (C.5).  Either  the  weighted  least 
squares  or  minimum  variance  solutions  may  then  be  obtained  by 
making  the  appropriate  substitution  for  the  matrix  M. 


bll  =  allmll  +  a21m21  +  a31m31 

b12  "  allm12  +  a21m22  +  a31m32 

b13  "  allm13  +  a21m23  +  a31m33 

b21  "  a12mll  +  a22m12  +  a32m31 

b22  *  a12m12  +  a22m22  +  a32m32 

b23  =  a12m13  +  a22m23  +  a32m33 
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G  * 


"-1 

C11 

bll 

b12 

b13 

gll 

g12 

g13 

-1 

f21 

C22^ 

b21 

b22 

b23_ 

_g21 

g22 

923_ 

gll  *  C11  bll  +  c12  b21 
g12  "  C11  b12  +  c12  b22 


g13 

g21 

g22 

g23 


C11  b13  +  c12  b23 
C21  bll  +  C22  b21 
C21  b12  +  C22  b22 
C21  b13  +  C22  b23 


Finally  we  obtain 


(C.  7) 


Ax  9h(tiq~t1p)  +  g12(T2Q_T2P)  +  g13(x3Q_T3P) 


LAy  921(t1Q_T1P)  +  922(t2Q"T2P)  +  g23 ( T3Q_T3P)  J 
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Appendix  D 


Calculation  of  R 


for  3  TD  Equal  Gaussian  Noise  Case 


From  Appendices  B  and  C,  R  was  defined  as 
(D.l)  R  =  <eeT> 

where  e  was  the  vector  of  time  difference  measurement  noises,  as¬ 
sumed  to  be  zero  mean.  Therefore, 


(D.  2) 


where  is  the  noise  on  the  ifc^  time  difference  measurement.  Now, 
since  each  time  difference  measurement,  T^,  is  essentially  the  re¬ 
ceipt  time  of  the  i^*1  slave  pulse  minus  the  receipt  time  of  the 
master  pulse,  then 


(D.  3) 


e 


i 


ni 


n 


m 


where  n^  is  the  tracking  filter  noise  jitter  on  the  receipt  time  of 
the  ifck  slave  pulse  and  nm  is  the  tracking  filter  noise  jitter  on 
the  receipt  time  of  the  master  pulse. 

Now  in  actuality,  each  has  additional  noise  added  to  it 
due  to  the  quantization  process  of  each  time  difference  counter.  It 
is  assumed  that  it  is  negligible.  There  may  even  be  additional 
smoothing  (averaging)  of  each  T\  after  the  counting  process.  In  any 
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event,  we  assume  that  the  relative  noise  performance  between 
Tl'  T2'  T3  rema*ns  the  same  through  these  additional  processes. 
Therefore,  it  is  sufficient,  for  G  matrix  purposes,  to  consider 
the  tracking  filter  noise  jitter  since  only  the  relative  noise 
statistics  will  affect  the  G  matrix. 


Assuming  <n^>  =  0,  <nm>  =  0,  then 


(D.  4) 


R=<ee  >  = 


R11  R12  R13 


R12  R22  R23 


R13  R23  R33 


where 

(D.  5) 


Rii  "  <Ei2>  *  <nt2>  +  <nm2>  1  “  1'2'3 


R.  .  =  <e  .  e  .  >  =  <n  > 
xj  x  ]  m 


i^j  i  -  1,2 
j  =  2,3 


since  it  is  assumed  that  <n^n^>  =  0.  Now  there  are  many  formulae 
for  calculating  noise  jitter  performance  out  of  a  tracking  filter. 
They  differ  according  to  the  order  and  bandwidth  of  the  tracking 
loop  and  the  specific  mechanization  of  the  phase  error  detector. 
We  assume  that  all  tracking  filters  are  identical  so  that 


(D.  6) 


n . 

l 


=  kQ/(S/N)i 


i  =  l,2,3,m 


where  kQ  is  a  function  of  the  specific  tracking  filter  implementa¬ 
tion  and  (S/N) ^  is  the  power  signal-to-noise  ratio  at  the  input 
to  each  tracking  loop.  Since  atmospheric  (external)  noise  is  dom¬ 
inant  (receiver  noise  negligible)  and  since  we  assume  the  RF  re¬ 
ceiver  input  bandwidths  are  all  the  same,  then  the  noise  presented 
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to  each  tracker  is  the  same  (even  if  the  RF  gains  from  antenna 
input  to  tracking  loop  input  is  different  between  receivers,  the 
signal-to-noise  ratio  remains  the  same) .  Therefore, 

_  k  N 

(D.7)  <n/>  *  i  =  1,2, 3,m 


so  that  the  relative  tracking  jitter  performance  is  only  a  func¬ 
tion  of  the  signal  powers,  S^,  into  the  receiver  antenna  input. 


(D.8) 


S. 

i 


P. 

i 


i  =  1,2,3 


and 


(D.  9) 


S 


m 


P 

m 


( 


4tt  r, 


2 

) 


where  P.^  is  the  effective  transmitter  power  of  the  ith  station; 
X=  wavelength  of  carrier  frequency;  r2i  and  ^  are  the  slant 
ranges  defined  in  Appendices  B  and  C.  Therefore, 


(D. 10) 


<n . 

l 


>  s 


k,  r_  . 
1  2i 


2  kl  ri 
<nm  >  =  ~  P 


P~ 

l 

2 


m 


i  =  1,2,3 


so  that 

(D. 11) 


i  =  1,2,3 
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The  value  of  k^  is  unimportant  except  that  it  would  be  good  prac¬ 
tice  to  set  it  so  that  the  elements  of  R  are  near  unity.  We  let 


(D. 12) 


r  r  -1 

R11 


.-1 


-1 


'12 


-1 


12 


-1 


22 


n  -1  „  -1 

R13  R2  3 


-n 


13 


-1 


'23 


-1 


33 


then 

(D. 13) 


=  R11  R22  R33  +  2  R12  R23  R13 


“ 

R13 

‘  R22 

R12  R33 

Rn_1 

= 

(R22 

R33 

-r23  ) /A 

R12_1 

= 

(R13 

R23 

_R33R12^ 

R  -1 
R13 

= 

(R12 

R23 

-R22Ri3)/A 

R  -1 
R22 

= 

(RU 

R33 

-R132)/A 

R  -1 
R23 

= 

(R12 

R13 

-rhr23)/a 

R  -1 

L  33 

= 

(RU 

R22 

-R122)/a 

Note  that  since  the  elements  of  R  are  directly  proportional  to  k. , 
then  the  elements  of  R  are  inversely  proportional  to  k^.  For 
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this  case,  from  Appendix  A, 

„  ,  .T  - 1.  .  -1  T  -1 

G  =  (A  R  A)  AR 

T  -1  T  -1 

Therefore,  A  R  is  inversely  proportaional  to  k^  while  (A  R  A) 

is  directly  proportional  to  k^f  resulting  in  G  being  independent 

of  k^. 
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TO:  J.  G.  Wall 

FROM:  L.  J.  Levy  and  V.  Schwab 

SUBJECT:  Methods  for  Obtaining  Approximate  and  Exact  Solutions 

to  a  LORAN  Navigation  Equation 

REFERENCE:  [1]  APL/JHU  S2R-79-008,  G-Matrix  Analysis  for  USCG 

LORAN-C  HHE  Navigation  Program,  by  L.  J.  Levy  and 
V.  Schwab,  January  19,  1979 

The  G-Matrix  technique  proposed  by  the  U.S.C.G.  is  a 
computational  method  for  converting  two  or  more  measured  LORAN 
TDs  into  horizontal  coordinates  of  navigator  position.  It  yields 
a  linear  estimate  of  position  referenced  to  the  surveyed  values  of 
a  'near'  waypoint  or  track  point  and  employs  a  precomputed  G-matrix 
to  transform  the  TDs. 

In  the  winter  and  spring  of  1979  the  proposed  linear 
G-Matrix  scheme  as  applied  to  the  LORAN  mini-chain  in  the  St  Marys 
River  was  simulated  on  the  IBM  3033  computer  at  APL.  It  was  then 
discovered  that  the  navigation  errors  resulting  from  the  high  grid 
curvatures  were  much  too  large  to  be  tolerated;  a  modification  of 
the  linear  G-Matrix  technique  was  then  sought  which  would  reduce 
the  navigation  error  to  an  acceptable  level. 

Two  modifications  were  studied.  The  first  of  these  involved 
the  iterated  application  of  a  second-order  correction  to  the  linear 
G-matrix  solution.  The  second-order  correction  greatly  reduced  the 
first-order  navigation  error,  but  the  residual  errors  were  still 
above  an  acceptable  level.  It  was  then  apparent  that  the  exact 
solution  of  the  LORAN  TD  equation  was  needed  and  also  within  easy 
reach  if  an  iterative  process  for  obtaining  it  were  admitted.  The 
method  for  obtaining  the  exact  solution  which  was  devised  retains 
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the  linear  G-Matrix  solution  for  its  starting  estimate  of  navigation 
position.  It  therefore  also  retains  the  important  feature  of  the 
surveyed  TDs  and  locations  of  the  waypoints. 

The  remaining  five  sections  of  this  memorandum  elaborate 
on  the  remarks  above.  Section  1  contains  a  derivation  of  the  basic 
LORAN  equation  relating  the  TD  vector  and  navigator's  position  vector. 
Section  2  defines  the  linear  G-Matrix  solution  of  the  LORAN  equation. 
Section  3  contains  a  derivation  of  the  second-order  correction  to 
the  linear  G-Matrix  solution.  The  exact  solution  is  presented  in 
Section  4.  The  three  methods  of  solution  are  compared  quantita¬ 
tively  in  Section  5  where  they  are  applied  to  navigation  in  the 
LORAN-C  mini-chain  of  the  St.  Marys  River. 
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1.  General  Form  of  the  LORAN  Time  Difference  Vector 

In  this  section  we  shall  develop  the  general  form  for  the 
LORAN  time  difference  vector  upon  which  the  three  navigation  solu¬ 
tions  derived  in  subsequent  sections  are  based. 

The  cartesian  coordinates  of  the  LORAN  receiver  aboard  a 
ship  are  denoted  by  X  and  Y  and  these  are  the  components  of  the 
vector 


z 


(1.1) 


Either  two  or  three  time  differences  (TDs)  are  measured 
and  for  the  sake  of  being  definite  we  shall  assume  in  this  memoran¬ 
dum  the  number  is  three.  The  time  difference  or  TD  vector  is  given 
by 


TD 


TD, 


TD. 


TD. 


(1.2) 


The  slant  ranges  from  the  ship  to  the  LORAN  slave  and 
master  stations  are  given  by 


[(X-X^2  +  (Y-Yj)2)*,  i=l,2,3,m 


(1.3) 
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in  which  isl,2 , 3  denotes  the  three  slaves  and  i=m  the  master  LORAN 
station. 

The  part  of  the  time  difference  vector  due  to  geometry 
alone  (i.e.,  the  ideal  part,  ignoring  biases  and  noise)  is  denoted 
by  h(je)  and  given  by 


<rrV/v 

hl 

h(z)  * 

tr2-rm)/v 

h2 

_(r3-rm)/v_ 

_h3_ 

(1.4) 


in  which  v  denotes  the  propagation  speed  of  the  LORAN  signal. 

The  matrix  of  partial  derivatives  of  the  elements  of  h 
with  respect  to  X  and  Y  is  denoted  by  A  and  given  by 


3hx 

3hx 

^x" 

W~ 

3h2 

3h2 

3X- 

■3Y" 

3h3 

3h3 

’3X“ 

W 

Let  2^  and  z^  represent  two  different  values  of  z.  The 
vector  difference,  h(Zg)  -  h(Zp)  ,  can  be  resolved  into  a  linear  com¬ 
ponent  equal  to  the  vector  A^tZg  -  and  a  residual  non-linear 
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component ,  denoted  by  the  vector  f  (Zp ,  z^) .  The  non-linear  component 
will  be  defined  by 


“  2lt Zq)  -  h (Zp)  -  AptzQ-Zp] 


(1.6) 


At  a  waypoint  P  a  measured  TO  vector  is  represented  by 

TDp  *  h(Zp)  +  Wp  +  rip  (1.7) 

in  which  Wp  represents  biases  and  warpage  in  the  measurement  and 
rip  represents  zero-mean  noise.  Similarly,  at  an  aribtrary  nearby 
point  Q  the  measured  TO  vector  is  represented  by 

TPq  -  h (Zg)  +  Wq  +  2q.  (1»8) 


The  surveyed  TO  vector  at  z»z^  is  the  mean  of  many  TD 
measurements  of  the  form  (1.7)  and  is  denoted  by  TE_ .  Since  np  has 
zero  mean  we  have 


TDp  -  h(Zp)  +  Wp. 


(1.9) 


If  P  and  Q  are  sufficiently  close  it  may  be  reasonable  to  assume 
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If  we  assume  (1.10)  we  find  from  (1.8)  and  1.9) 


TDq  -  TDp  «  h(z^)  -  h(2p)  +  rig.  (1.11) 

However  by  means  of  (1.6)  we  are  able  to  eliminate  h(z) 
from  (1.11)  and  obtain 

TDq  -  TDr  -  Apl^  -  Zpl  +  f (Zp,^)  +  (1.12) 

In  (1.11)  and  (1.12)  Zp  represents  the  surveyed  position  of 
point  P.  The  difference  between  the  mean  of  many  TD  measurements 
made  at  P,  TDp,  and  h(zp)  is,  by  definition,  the  warpage  or  bias  at 
P  which  is  denoted  by  Wp  in  (1.7).  All  three  of  these  vectors  are 
known  for  each  waypoint. 

The  problem  of  LORAN  navigation  is  to  determine  the  position 
vector  Zq  given  the  measurement  vector  TDq.  The  different  ways  in 
which  this  has  been  done  are  explained  in  the  following  three  sections 


2.  Linear  G-Matrix  Navigation  Solution 

The  linear  G-Matrix  solution  to  the  LORAN  navigation  pro¬ 
blem  represented  by  (1.12)  is  obtained  under  the  assumption  that 
or 


2£q  -  lip  -  '+  V 


(2.1) 
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In  the  notation  of  Ref.  [1]  let 


u  ■  TD_  -  TD, 


X  -  Z-  -  2_ 
“0  HP 


e  =  *Q 


R  *  Cov(e) 


W  positive  definite,  symmetric  weighting  matrix. 


Then  (2.1)  takes  the  form 


u  «  Ax  +  e. 


(2.2) 


A  solution  of  (2.2)  is  represented  by 


in  which 


x  *  G  u 


G  «  (ATWA)_1ATW. 


(2.3) 


(2.4) 


Note  from  (2.4)  that  GA-X.  The  solution  given  by  (2.3) 
is  called  the  ordinary  least  squares  solution  when  W«I ;  the  weighted 
least  squares  solution  when 
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the  minimum  variance  solution  when  W=R  .  In  any  case  the  solution 

*  A 

of  (2.1)  is  denoted  by  z^  =  x  +  z^  and  is  given  by 

Iq  -  Zp  +  Gp  ITD Q  -  TDp] .  (2.5) 

3.  G-Matrix  Solution  with  2nd  Order  Correction 

In  Appendix  A  we  obtain  the  second-order  expansion  of  a 
vector  function  (TD)  of  two  variables  (X,Y) .  The  result  is  expressed 
by  (A.6).  When  we  apply  this  result  to  (1.12)  we  find 

ISq-ISp  “  W5*1  +  2Q 


(3.1) 


in  which 
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(3.3) 


and  Ap  is  defined  by  (1.5). 

The  solution  of  (3.1)  is  obtained  through  an  iterative 
process.  In  the  first  iteration  we  compute  the  linear  approxima¬ 
tion  defined  by  (2.5),  i.e.. 
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In  general. 


~n-l 

o 


-  (i/2)  i  *  [g 

i=l  -1  ^ 


n-1  ,T 

-Sp] 


BPi^Q  “-P1 


(3.8) 


and 


(3.9) 


The  iterative  process  may  be  terminated  when  the  change  in  a  from 
the  previous  iteration  becomes  negligible. 


4.  Exact  Solution  of  the  Non-Linear  LORAN  Equation 

We  shall  present  in  this  section  an  iterative  scheme  for 
obtaining  an  exact  solution  to  the  non-linear  LORAN  equation  given 
by  (1.12)  or 

TDq  -  TDp  «  Apl^-Zp]  +  f (ZqeZj,)  +  2q.  (4.1) 

Again  we  use  the  linear  G-matrix  solution  as  the  starting  estimate 
in  the  first  iteration.  Thus  we  have 

-  ip  +  S  'I2q  - 
Then  from  (1.6)  we  obtain 

I1  -  h(j£)  -  h(Zp)  -  Ap [ -  Zp).  (4.3) 
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For  the  second  Iteration  we  find 


z  +  G  [TD 
^P  P 


and 

i2  -  h<i|)  -  M*p>  -  *p[«|  -  if)- 


In  general. 


fn_1  *  h(z£_1) 


-  h(Zp)  - 


*P 


tiS"1 


-  *,) 


and 


A  simpler  scheme  than  that  represented  by  (4.6) 
can  be  obtained  by  avoiding  £  and  dealing  with  h  alone. 

A  1 

In  (4.4)  if  we  substitute  from  (4.3)  for  f  we 


3 


2l 


r^l 


£  '  Gp  [h(«J)  -  h(«p)]  +  GpAp  [«2  -  «p] 


or,  since  GpAp“ 


£q  “  21q  “  Ip  "  Gp  ^(Iq)  “  • 


(4.4) 

C4. 5) 

(4.6) 

(4.7) 

and  (4.7) 

obtain 


(4.8) 
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Por  the  third  iteration  we  find  from  (4.7) 


4  -  4  -  gp£  • 


(4.9) 


From  (4.9)  and  (4.5)  we  obtain 


ia  +  30  -  ip  -  S  -  il'ip”- 


(4.10) 


If  we  add  the  term 

A 1  A  1  _  .  a1  .  .  a!  . 

zj  -  +  Gp  h(z£)  -  Gph(zJ)  -  0 

to  the  right-hand  side  of  (4.10)  and  apply  (4.8)  to  the  result  we 
find 


or 


i! 


-  ip  -  '  + 

A O  A 1  a9  «  A 1  „ 

24  -  4  -  s^(4}  -  ^(4)]* 


(4.11) 
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In  general. 


zn 


-  %~2  -  -  h(z£“2)  ] 


(4.12) 


The  iterative  scheme  represented  by  (4.12)  may  be  terminated  when 

A 

the  change  in  z  from  the  previous  iteration  is  sufficiently  small 
to  be  negligible. 


5.  Comparison  of  the  Navigation  Solutions  Obtained  by  Means  of 

the  Linear,  Linear  witn  2ncl  Order  Correction  and  Exact  Methods 
of  Solving  the  LORAN  TD  Equation 

A  comparison  of  the  results  obtained  by  the  three  methods 
of  solving  the  LORAN  TD  equation  (1.12)  has  been  made  for  navigation 
between  selected  waypoints  in  the  LORAN-C  mini-chain  in  the  St.  Marys 
River  connecting  Lakes  Superior  and  Huron. 

In  Table  5-1  the  integers  NWAY  denote  waypoint  numbers;  S 
represents  the  distance  from  a  waypoint  P  to  a  point  Q  on  the 
straight  line  joining  waypoint  P  with  the  next  waypoint.  When  S=0, 

Q  is  at  P;  when  S  attains  its  maximum  value  at  a  waypoint  (e.g., 
S=11.19  at  waypoint  1)  Q  is  at  the  next  waypoint.  ATE  and  CTE  denote 
the  along-track  error  and  cross-track  error  in  navigated  position  at 
point  Q. 

Each  row  in  Table  5-1  represents  a  navigation  run.  The 
distance  S,  of  point  Q  from  Point  P  is  increased  until  the  next 
waypoint  is  reached. 

The  three  double  columns  of  ATE  and  CTE  contain  the  naviga¬ 
tion  errors  which  result  using  the  linear  G-Matrix  solution  defined 
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in  Section  2;  the  iterative  G-Matrix  solution  with  second-order 
correction  defined  in  Section  3 ;  the  iterative  exact  solution  de¬ 
fined  in  Section  4. 

It  may  be  seen  from  Table  5-1  that  the  navigation  errors 
increase  with  increasing  S  except  in  the  case  of  the  exact  method 
of  solution.  Also,  while  the  second-order  correction  achieves  a 
considerable  reduction  in  navigation  error  over  the  straight  linear 
method,  the  errors  remain  unacceptably  large  especially  when  the 
distance  between  waypoints  is  great. 
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TABLE  5-1. 

COMPARISON  OF 

METHODS  OF 

SOLUTION 

Exact 

Linear 

2nd  Order 

Correction 

S  (KYDS) 

ATE  (FT) 

CTE  (FT1 

ATE  (FT) 

CfE  (rTJ 

ATE  (FT) 

CTF  (FI) 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

O’.O 

1.1C0CU 

-1 19.0 

10.0 

0.6 

-0.1 

0.05 

-0.01 

2.20000 

-468.6 

35.2 

5.1 

-1.3 

0.03 

-0.00 

3.30CC0 

-1044.3 

86.1 

19.0 

-4.5 

-0.11 

0.01 

4  .40000 

-1824.7 

148.5 

48.5 

-11.3 

0.28 

-0.03 

5.50000 

-2627.4 

225.4 

ICO. 2 

-22.9 

C.21 

-0.02 

6. 600C0 

-4008. 1 

313.1 

180.2 

-40.6 

0.20 

-0.02 

7.70000 

-5361.6 

405.6 

254.6 

—66.3 

0-23 

-0.02 

e.eocoo 

-6871.5 

512.1 

448.5 

-ICO. 7 

0-26 

-0.02 

9 .50000 

-8522.5 

617.  5 

646.5 

-145.0 

0.37 

-0.03 

11.000C0 

-10298.7 

724.5 

892.0 

-200. 9 

0.4$ 

-o.ou 

11.  1S0  24 

-10617.6 

742.8 

935.6 

-211.2 

-0.37 

0.03 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

i. ecoco 

-48.2 

-25.5 

-0.5 

-0.8 

0.0  0 

0.00 

2. 000  CO 

-155.0 

-105.2 

-3.7 

-6.3 

-C.01 

-0.01 

3. coooo 

-443.4 

-244.2 

-11.  7 

-20.7  | 

0.02 

0.01 

4.  ocooo 

-756.  1 

-447.4 

-26.0 

-4  8.1 

-0.03 

-0.02 

5.00000 

-1255.9 

-715.4 

-47.6 

-91  .8 

C-OU 

0.03 

5. 417  SI 

-1480.3 

-854.6 

-58.9 

-113.6 

0.08 

0.05 

c.o 

0.0 

O.C 

0.0 

0.0 

-O-CO 

-0.00 

C.  863  79 

-20.9 

-23.2 

-0.0 

-0.6 

-0.02 

-0.02 

0.0 

0.0 

C.O 

0.0 

0.0 

-0.00 

-0.00 

1. COOOO 

11.5 

7.8 

0.4 

-0.4 

0.00 

0.00 

1.56471 

28.6 

ie.6 

1.7 

-1.5 

0.01 

0.02 

0.0 

0.0 

0.0 

0.0 

0.0 

-0.00 

-0.00 

1.00000 

-11.5 

-31.6 

0.2 

-1.1 

0.00 

0.00 

2.00000 

-45.1 

-131.0 

2.0 

-8.9 

C.  03 

O.OU 

3.00CC0 

-55.5 

-305.2 

7.1 

-25.7 

-0.C2 

-0-03 

4.000  CO 

-173.5 

-560.5 

18.0 

-69.3 

0.02 

0.02 

4.2E515 

-198.0 

-649.6 

22.5 

-64.6 

0.  C  3 

0.C4 

0.0 

0.0 

0.0 

0.0 

0.0 

-0.00 

-0.00 

1. COOCO 

-39.2 

-35.1 

0.4 

0.4 

0.00 

0.00 

2 .COOOO 

-135.1 

-138.7 

3.1 

3.3 

0.  06 

0.02 

3. 000  CO 

-344.6 

-307.0 

10.6 

12.0 

-0.02 

-0.01 

3. 55559 

-450.5 

-436.4 

10. S 

21.6 

-O-Or 

-0.04 

0.0 

0.0 

0.0 

0.0 

3.0 

0.0 

0.0 

1. OOUCC 

-34.5 

-12.4 

0.3 

1.0 

-0.02 

-o.co 

2.00000 

-  125.2 

-46.0 

5.9 

7.  7 

0.01 

0.  00 

3.00000 

-297.4 

-9  5. 4 

19.3 

2  5.6 

-C.00 

-0.00 

2. 29662 

-336. 7 

-112.4 

25.4 

33.9 

-O.C 

-0. 00 

0.0 

0.0 

0.0 

0.  0 

0.0 

0.0 

0.0 

C.  50921 

-4.  3 

13.3 

0.2 

0.2 

-0-00 

0.00 

0.0 

0.0 

0.0 

1.0 

3.0 

0.0 

0.0 

1. CCOCO 

-12.3 

16.6 

0.3 

0 . 4 

-C.00 

0  •  00 

2.000 CO 

-48.1 

67.6 

2.5 

2.o 

0.00 

-0.00 

2.C26C7 

-49.4 

69.4 

2.6 

2.5 

0.00 

-  j .  00 

ifc* 
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APPENDIX  A 

Expansion  to  the  2nd  Order  of  a 
Vector  Function  of  Two  Variables 


To  the  second  order  in  h  and  k  the  Taylor  expansion  of 
the  scalar  function  f(x,y),  where  x  and  y  are  independent  variables, 
is  given  by 


f (z+h,y+k)  =  f (x,y)  +  (hfx(x,y)  +  k  fy(x,y)) 

+  (1/2!)  jh2fxx(x,y)  +  2hkfxy(x,y)  +  k2fyy(x,y)J. 

(A.  1) 

We  define  the  1x2  matrix  A,  the  2x2  matrix  B  and  the  2x1  vector  v 
by  the  expressions  in 


A  = 


—  ** 

f  f 

r—  — 

xx  xy 

V  = 

h 

Hi 

Hi 

k 

L  Ky  yy  - 

(A.  2) 


The  matrices  and  vector  defined  in  (A. 2)  permit  us  to  write  (A. 1) 
in  the  form 


f (x+h,y+k)  ■  f(x,y)  ■  Av  +  (1/2!)  vTBv 


(A.  3) 
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To  verify  (A.  3)  observe  from  (A.  2)  that 


AV  -[hfx  +  kt  3 


BV 


h£xx  +  kfxy 


l!*£xy  +  kfyyj 


vTBv  «  [h  k) 


hf 

+ 

kf 

XX 

xy 

hf 

♦ 

X" 

Ml 

xy 

yy 

^'xx  +  2hkfxy  +  k2fyy 


and  note  that  these  results  agree  with  (A.  1) . 

Consider  next  the  3x1  vector  function  F (x,y)  defined  by 


F(x,y) 


fjtxsy) 

f2<x»yJ 

f3<x,y) 


(A.  4) 


Again  we  wish  to  form  the  second-order  Taylor  expansion  of  F(x,y). 
For  this  purpose  we  define  the  3x2  matrix  A,  the  2x2  matrices  B^, 
1*1, 2, 3  and  the  3x1  vector  i*l,2,3  in 
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A  -  f 


fixx  f ixy 


fixy  fiyy 


,  1-1,2, 3 


(A.  5) 


Then,  to  the  second  order  in  h  and  k. 


F (x+h,y+k)  «=  F(x,y)  ■  Av  +  (1/2J)  J  ti.vTB.v  (A. 6) 

i-1 


where  v  has  been  defined  in  (A. 2) . 


C-A-3 


TM(  JOHNS  HOSKINS  UNIVfNSITV 

APPLIED  PHYSICS  LABORATORY 

LAijaci.  ujavtJNO 


S2R-79-166 

ZT70S2R0 

October  22,  1979 


TO: 


J.  G.  Wall 


FROM:  L.  J.  Levy/V.  Schwab 

SUBJECT:  Design,  Simulation  and  Performance  of  Digital  TD  and 

Heading  Filters  for  the  USCG  LORAN-C  HHE  Navigation 
Program 

REFERENCES:  [1]  APL/JHU  TYA-0-002,  Compilation  of  Characteristics 

of  the  Second-Order  Digital  Filter,  by  R.  J. 
McConahy,  August  29,  1963 

[2]  IRE  Transactions  PGAC,  Vol.  AC-7,  July  1962, 
Synthesis  of  an  Optimal  Set  of  Radar  Track-while- 
scan  Smoothing  Equations,  by  T.  R.  Benedict  and 
G.  W.  Bordner 

[3]  APL/JHU  S2R-79-162,  Methods  for  Obtaining  Approxi¬ 
mate  and  Exact  Solutions  to  a  LORAN  Navigation 
Equation,  by  L.  J.  Levy  and  V.  Schwab,  October  8, 


The  time  difference  (TD)  measurements  made  in  the  USCG 
LORAN-C  HHE  Navigation  Program  may  be  expected  to  contain  zero-mean 
noise  errors  ranging  in  size  up  to  tens  of  nanoseconds.  Since  the 
TD  measurements  are  converted  directly  to  display  ship's  position 
on  a  scope,  it  is  desirable  to  filter  out  large  fluctuations  in  the 
TDs  before  computing  position. 

In  the  absence  of  ship's  maneuver  (turning)  the  TDs  can  be 
effectively  filtered  by  a  conventional  second-order  digital  filter 
which  has  zero  velocity  lag.  However  large  velocity  lags  and  posi¬ 
tion  errors  are  introduced  by  this  filter  during  the  frequent  turns 
required  on  the  St.  Mary's  River.  It  has  been  demonstrated  that  this 
deficiency  of  the  conventional  filter  can  be  essentially  eliminated 
by  the  addition  of  a  first-order  filter  which  provides  a  filtered 
value  of  acceleration  normal  to  the  ship's  course.  The  raw  accelera¬ 
tion  measurement  is  obtained  from  differential  measurements  of  ship's 
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heading  angle  which,  in  turn,  measure  the  turning  rate  of  the  ship. 

The  main  body  of  this  memorandum  is  divided  into  three 

parts : 

I  Filter  Design 

II  Filter  Simulation 

III  Filter  Performance 

The  section  on  filter  design  is  divided  into  two  parts. 

The  first  part  reviews  the  theory  of  the  conventional  second-order 
digital  a,b-filter;  the  modification  to  the  filter  required  when 
aged  measurements  and  display  delay  are  present;  the  reduction  of 
velocity  lag  by  means  of  acceleration-aided  prediction  in  the  filter. 
The  general  filter  equations  developed  in  the  first  part  are  applied 
in  the  second  part  to  the  particular  problem  of  filtering  LORAN  TD 
measurements . 

The  section  on  filter  simulation  contains  a  fairly  de¬ 
tailed  description  of  a  digital  simulation  prepared  for  the  IBM  3033 
Computer  at  APL  in  the  spring  of  1979.  The  purpose  of  the  simulation 
was  the  provision  of  the  possibility  of  the  computation  of  filter 
performance  in  order  to  confirm  the  choices  of  filter  gains  and 
other  optional  system  parameters. 

The  section  on  filter  performance  presents  a  small  repre¬ 
sentative  sample  of  the  simulation  runs  which  were  made  on  the  3033 
Computer.  The  particular  series  of  four  runs  were  chosen  to  show 
the  effect  on  filter  performance  of  variations  in  the  gain  of  the 
heading  (acceleration)  filter.  These  runs  show  clearly:  (1)  the  need 
to  greatly  reduce  velocity  lag;  (2)  the  ability  of  the  heading  filter 
to  essentially  eliminate  velocity  lag. 


LJL/VS : nt 
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1-1.1  The  Conventional  a,b-Filter 

The  following  description  of  a  second-order  digital  filter 
used  for  filtering  and  estimating  the  rate  of  change  of  a  digital 
data  sequence  follows  closely  the  excellent  account  given  in  Ref  Ill. 
The  sequel  is  often  nearly  but  not  precisely  verbatim  and  the 
notation  of  the  present  text  also  differs  slightly  from  that  of 
Ref  HI. 

Suppose  one  has  a  sequence 
x* (n) ,  n*0 ,1,2,. . . 

of  input  or  raw  data.  We  suppose  this  sequence  is  measured  (or  de¬ 
livered)  at  a  sample  period  of  At  seconds.  Then  x*(n)  would  be  the 
measurement  of  the  quantity  'x'  at  time  t=nAt. 

Obviously  the  individual  measurements  will  contain  noise. 

It  is  often  desirable  to  smooth  them  so  that  the  system  functions 
which  use  the  data  will  not  operate  wildly  and  erratically,  as  they 
might  if  nothing  were  done.  But,  more  than  this,  there  is  usually 
a  desire  to  estimate  the  rate  of  change,  x,  of  the  quantity  x.  This 
is  very  useful  in  case  there  is  a  sudden  fade  in  the  data  sequence 
x*(n)  and  we  are  forced  to  ’coast'  awhile. 

One  of  the  simplest  filtering  or  estimation  schemes  that 
will  allow  us  simultaneously  to  smooth  the  data  and  estimate  rate  of 
change  is  the  second-order  filter  defined  by 


r(n) 

-  x* (n) 

i 

x  > 

:T 

x  > 

-  x'  (n) 

+  a  r(n) 

a) 

*(n) 

*  x' (n) 

+  !tr(n) 
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and 

^  a  a 

x'(n+l)  «  x(n)  +  x(n)At 

:  :  (2> 

x'  (n+1)  *  x(n) 

in  which 

x(n)  is  the  smoothed  estimate  of  x  at  step  n 

a 

x(n)  is  the  smoothed  estimate  of  x  at  step  n 

x' (n)  is  the  predicted  value  of  x  at  step  n  and  usually 
made  at  step  n-1 

A 

x'  (n)  is  the  predicted  value  of  x  at  step  n  and  usually 
made  at  step  n-1 

a  position  smoothing  gain 

b  rate  smoothing  gain 

The  equations  in  (1)  update  the  predictions  (primed  quanti¬ 
ties)  made  for  step  n  to  estimates  of  x  and  x.  The  quantity  r(n)  is 
the  error  in  the  prediction  of  x(n)  referenced  to  x*(n).  Note  that 
if  a«0,  x(n)»x'(n)  and  no  weight  is  given  to  the  measurement  x*(n); 
however  if  a*l,  x(n)«x*(n)  and  the  measurement  at  step  n  is  given 
full  weight.  The  equations  in  (2)  define  the  prediction  for  one  time 
step  in  advance. 

Equations  (1)  and  (2)  can  be  reduced  to  the  second-order 
difference  equation 

x(n+2) -(2-a-b)  x(n+l)  +  (l-a)x(n)  *  ax*(n+2)  +  (b-a)  x*  (n+1) . 

The  performance  of  the  second-order  filter  depends  on  the  choice  of 
the  smoothing  gains  a  and  b. 

A 

For  stability  of  x(n)  and  x(n)  it  is  required  that 
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b  <  4  -  2a,  0  <  a  <  2. 

The  filter  is: 

2 

underdamped  if  4b  >  (a+b)  ; 

2 

critically  damped  if  4b  *  (a+b)  ; 
overdamped  if  4b  <  (a+b)  . 

The  filter  is  said  to  be  optimum  if  a  and  b  satisfy  the 

relation 


which  results  in  a  slightly  underdamped  filter. 

We  shall  now  explain  what  is  meant  above  by  the  term  'opti¬ 
mum  filter'.  Assume  the  noise  on  each  input  sample  has  zero  mean, 

2 

variance  (mean  square)  oN  and  that  the  noise  on  different  samples  is 
uncorrelated.  The  ability  of  the  filter  to  smooth  the  data  is  fre¬ 
quently  measured  by  the  index 

r  (o)  «  Variance  in  smoothed  position  output 
x v  '  Variance  in  raw  position  input 

which  is  referred  to  as  the  variance  reduction  ratio.  There  is  also 
a  variance  reduction  ratio  for  the  smoothed  output  velocity  sequence, 
Jt(n)  ,  which  is  denoted  by  R^(0).  For  any  stable  combination  of  a  and  b 


v°> 


2a2+  b (2  -  3a) 
a  ( 4  -  b  -  2a) 
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Ri(0> 


It  is  a  general  characteristic  of  the  second-order  filter 
that,  given  sufficient  time,  it  will  settle  onto  a  linear  input 
sequence,  i.e., 


x*(n)  *  An  +  B 


with  negligible  error.  That  is 


lim  (x(n)  -  x*(n) ]  =  0 


for  such  a  linear  sequence.  This  property  is  usually  referred  to 
as  'zero  velocity  lag'. 

A  useful  index  is  one  based  on  the  unit-ramp  response.  L 


h^Cn)  *  x(n)  for  x*(n)  =  n 


assuming  the  filter  is  initially  at  rest.  Then  the  index  is 


C  -  l  (n  -  h  (n) )  . 


One  can,  similarly,  define  a  velocity  index 


C-  -  l  (1  -  h*  (n) ) 
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where 

h^n)  *  x(n)  if  x*(n)  *  n. 

For  any  stable  combination  of  a  and  b 

C  =  (2-a)  (1-a)  2 
x  ab(4  -  b  -  2a) 

1  .a2  (2-a)  +  2b  (1-a) 

,  a  - j  -  . 

x  At  ab  ( 4  -  b  -  2a) 


The  indices  and  Rx(0)  each  depend  on  the  variables  a 
and  b,  and  one  might  compute  contours  of  these  functions  in  the  a,b- 
plane  using  the  formulas  given  above.  The  following  problem  in  opti¬ 
mization  is  posed  by  the  two  indices  and  may  be  formulated  in  terms 
of  two  questions: 


(i)  For  what  values  of  a  and  b  is  Cx  minimized  for  a  given 
desired  value  of  Rx(0)? 

(ii)  Conversely,  for  what  values  of  a  and  b  is  Rx(0)  minimized 
for  a  given  desired  value  of  Cx(0)? 

The  same  two  questions  can  be  asked  with  the  velocity  index 
R* (0)  replacing  R  (0)  and  C*  replacing  C  .  It  is  proved  in  Ref.  [2] 
that  all  four  questions  have  the  same  answer:  b**av(2-a)  >  as  asserted 
above.  It  is  in  this  sense  then  that  the  relationship  of  position 
and  rate  gains  stated  above  define  an  optimum  filter. 
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In  Table  1-1  there  are  listed  some  characteristics  first, 
of  a  critically  damped  and  then  an  optimum  second-order  filter.  In 
each  case  six  values  of  the  position  gain  are  specified  in  column  1 
and  the  corresponding  value  of  the  rate  gain,  b,  computed  by  means 
of  the  formula  given  at  the  top  of  the  table,  appears  in  column  2. 
The  quantities  0pos  and  0vEL  are  the  square  roots  of  Rx(0)  and  Rx(0) 
defined  above.  The  values  of  the  filter  time  constant  were  deter¬ 
mined  using  formulas  for  the  position  response  function  given  in 
Ref.  [1]  but  not  reproduced  here. 

1-1.2  Conventional  Filter  with  Aged  Measurement  and  Display  Delay 

Assume  that  at  time  tn=nAt  the  measurement  x*  has  an  age 
At.  >0  (i.e.,  the  measurement  is  made  at  t  -At-) .  Also  assume  that 
the  measurement  will  not  be  displayed  until  after  a  delay  At^M)  at 
tn+AtD*  The  equations  for  updating  the  estimate  at  tR-AtA  are 


r(n-AtA) 

x(n-Atft) 

x(n-At. ) 
A 


x*(n-AtA>  -  x’(n-AtA) 
x'(n-At A)  +  a  r(n-Atft) 
x'(n-AtA)  +  ^  r(n-AtA). 


(3) 

(4) 

(5) 


The  forward  predictions  to  the  next  update  time  tn+1-AtA  are  made 
by 


x'(n+l-AtA)  ■  x(n-AtA)  +  x(n-AtA>At 


(6) 


x  '(n+l-AtA) 


x(n-AtA> 


(7) 
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The  forward  predictions  to  the  display  time  tn+ Atp  are  made  by: 


xD (n+AtD) 

*  x(n-AtA)  +  x(n-AtA) (AtA+At^) 

(8) 

xD  (n+Atp) 

A 

=  x(n-AtA) . 

(9) 

We  shall  next  arrange  the  equations  above  so  that  updating 

is  done  at  t  +At_.  First  we  evaluate  (8)  and  (9)  at  t  ,,+At_  and 
n  D  n+ 1  d 

obtain 


xD(n+l+AtD)  =  x(n+l-AtA)  +  x(n+l-AtA) (AtA+AtD) 
xD(n+l+AtD)  -  x(n+l-AtA) . 

Next  we  evaluate  (3) ,  (4)  and  (5)  at  tR+1~AtA  and  find 

r(n+l-AtA)  =  x*(n+l-AtA)  -  x' (n+l-AtA) 

x(n+l-AtA)  .=  x' (n+l-AtA>  +  a  r(n+l-AtA) 

A  A  . 

x(n+l-AtA)  =  x'  (n+l-AtA)  +  (n+l-AtA) 


(10) 

(11) 


(12) 

(13) 

(14) 


By  means  of  (13)  and  (14) ,  (10)  and  (11)  become 

x°(n+l+At  )  =  x' (n+1-At.)  +  a  r(n+l-At  ) 

D  A  A 

+  [x'  (n+l-At. )  +  yr-r  (n+l-At.  )]  (At. +  Atn)  (15) 

A  At  A  A  D 

xD(n+l+At  )  =  x'  (n+l-At.)  +  7—r (n+l-At. ) .  (16) 

D  A  At  A 
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Then,  by  means  of  (6)  and  (7) ,  (12) ,  (15)  and  (16)  be¬ 
come 

r(n+l-dtA)  ■  x*(n+l-dtA)  -  x(n-AtA)  -  x(n-AtA)At  (17) 

xD(n+l+AtD)  -  x(n-AtA)  +  x(n-AtA)At  +  a  r(n+l-AtA) 

+[x(n-AtA)  +  |^r(n+l-AtA)] (AtA+AtD)  (18) 

xD(n+l+AtD)  »  x(n-AtA)  +  —r  (n+l-AtA) .  (19) 

Equations  (8)  and  (9)  can  be  used  to  express  (17) ,  (18)  and  (19)  in 
the  forms 

r(n+l-AtA)  -  x*(n+l-AtA)  -  xD(n+AtD)  +  xD(n+AtD) (AtA+AtD) 

-  xD(n+AtD)At 

A 

-  x*(n+l-AtA)  +  iD(n+AtD)  (At^Atp) 

-  (x^n+Atp)  +  iD(n+AtD)At] .  (20) 

A 

xD(n+l+AtD)  -  xD(n+AtD)  +  xD(n+AtD)At 

+  [a  +|^(AtA  +  AtD) ]  r(n+l-AtA)  (21) 

xD(n+l+AtD)  -  x(n+Atp)  +  ^r(n+l-AtA)  (22) 
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Let 


At 


AD 


AtA  +  AV 


(23) 


The  predicted  values  at  time  tn+^+AtD  are 


xD  (n+l+AtD)  «  xD(n+AtD)  +  xD(n+AtD)At 


xD  (n+l+Atp)  *  fcD(n+AtD) 


(24) 

(25) 


Then  by  means  of  (23),  (24)  and  (25)  equations  (20),  (21)  and  (22) 
become 


r(n+l-AtA)  *  x*(n+l-AtA)  +  fcD(n+l+AtD)  At^ 


-  xD  (n+l+AtD) 


XD(n+l+AtD)  «  xD  (n+l+^tp)  +  [a  +  |tAtADlr(n+1"AtA) 


(26) 

(27) 


kD(n+l+Atjj)  »  (n+l+AtD)+  ~r(n+l-AtD)  . 


(28) 


1-1.3  Addition  of  a  Filtered  Acceleration  Measurement  to  Aid 
Predi ctlon 

While  the  second-order  filter  defined  by  (26)  ,  (27)  and 
(28)  is  a  zero-velocity  lag  filter  it  will  develop  a  lag  in  the  pre¬ 
sence  of  a  non- zero  acceleration,  Sc.  If  acceleration  can  be  measured 
and  introduced  into  (26)  by  the  addition  of  a  term 


0.5  AtAE)X 
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the  lag  due  to  acceleration  can  be  greatly  reduced.  However  a 
noisy  measurement  of  acceleration,  x* (n) ,  would  first  have  to  be 
filtered.  To  accomplish  the  filtering  we  propose  the  first-order 
acceleration  filter  represented  by 

£<n)  »  £(n-l)+  c[x* (n)  -  X(n-D).  (29) 

Note  in  (29)  that  if  the  gain  c«0,  no  weight  is  given  to  the  accel- 

V\  A  A  . 

eration  measurement  and  x(n) «x(n-l) «0;  however,  if  c«l,  x(n)«x  (n) 
and  the  measurement  receives  full  weight. 

1-2.  Application  of  the  Acceleration- Aided  Second-Order  Filter  to 
tne  Problem  of  Filtering  LORAN  TD  Measurements 

In  Ref.  [3]  the  LORAN  TD  (time  difference)  is  given  in 
'exact'  form  by  Eq.  (1.12)  and  in  linearized  form  by  Eq.  (2.1)  or 
the  expression 


TDq  -  TDp  +  Ap [Zg  -  Zy]  +  ng 

in  which  TD  denotes  a  vector  of  two  or  three  TDs ;  Z_  denotes  a  two- 
dimensional  horizontal  position  vector  of  the  navigator;  n  is  a 
vector  of  zero-mean  noise  errors  in  the  TD  measurement;  A  is  a  matrix 
of  first-order  partial  derivatives  of  the  TDs  with  respect  to  the 
navigator's  position  coordinates.  Subscript  P  denotes  a  variable 
whose  value  is  known  and  associated  with  a  fixed  waypoint  with  co¬ 
ordinates  Z»Z_ .  Subscript  Q  denotes  a  variable  associated  with  the 
variable  point  Q  with  the  generally  unknown  coordinates  £”Zq*  From 
the  measurement  TDq,  a  'navigated'  value  Zq  is  to  be  inferred. 

If  we  consider  a  single  TD  rather  than  two  or  three  as  in 
the  vector  equation  above  we  can  use  the  equation 
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TDq  «  TDp  +  a^txQ-Xp)  +  *12(yQ-yp)  ♦  »  (30) 

in  which  all  quantities  are  scalars. 

There  are  four  time  intervals  which  appear  explicitly  in 
the  equations  representing  the  TD  and  heading  angle  (i.e.,  accelera¬ 
tion)  filters.  These  are: 

ATK  Time  interval  between  TD  measurements 

ATl  Time  interval  between  heading  measurements 

dT.  Incremental  time  step  in  dead  reckoning  computation 
1  in  TD  filter 

ATA  Time  delay  in  the  measurement  of  TDs. 

In  addition/  there  is  a  time  delay,  dT ^ ,  in  the  measurement  of  ship's 
heading  angle;  however  this  quantity  does  not  appear  explicitly  in 
the  filter  equations  which  follow. 

In  the  time  interval  between  TD  measurements,  i.e.,  when  a 
TD  measurement  is  not  available,  predictions  are  made  by  means  of  the 
formulas 

A  A  ^ 

TDi  =  TDi_1  +  TDi_1dTi  +  0.5  TD  dT^  (31) 

A  A  £ 

TDi  «  TD^  +  TD  dT^.  (32) 

When  a  TD  measurement,  denoted  by  TD^,  is  available  the 
following  computations  are  made: 


■  TDi 

(33) 

A 

■  fDi 

(34) 
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r,.  -  TD*  +  ATa  TDjJ  +  0 . 5 AT ^  TD  -  TDV’ 


(35) 


TDk  -  TDk  +  [a  +  b(ATA/ATK)lrk 
T*Dk  -  TDk  +  (b/ATK)rk. 

The  function  of  the  heading  filter  is  to  provide  the 
filtered  term  TD  appearing  in  (31)  and  (35) .  Prom  (30)  we  have 


(36) 

(37) 


TD  -  axl*  +  a12J 


and  by  geometry 

x  *  v  sin  ij) 

$  »  v  cos  tp 

where  v  denotes  the  constant  (i.e. ,  slowly  changing)  speed  of  the 
ship  and  ip  is  ship's  heading  angle  measured  clockwise  from  north. 
The  accelerations 

x  «  VTp  cosiJ> 
y  ■  v^  sintp 

can  be  approximated  using  measurements  ip*  and  ®ade  ATL  seconds 

apart.  The  normal  acceleration  measurement 

<  ■  v<*n  -  Cl1/AT1 
is  filtered  in  the  first-order  filter 
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Then 


^  ^  .  e  a 

an  ’  an-l  +  clan  '  an-l! • 

TD  -  |«ucos*J  -  «12,inO  V 


(38) 


(39) 


Equations  (31)  through  (39)  represent  the  TD  and  Heading 
Filter  applied  to  a  single  TD.  These  equations  are  based  on  the 
second-order  filter  with  acceleration  aiding  which  is  described  in 
Section  1-1.  The  equations  are  used  in  the  filter  simulation  program 
described  in  Section  II. 
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Table  1-1 

Summary  of  Some  Characteristics  of 
The  Second-Order  Filter 


Critically  Damped  Filter:  4b  ■  (a+b)2 


a 

b 

°POS 

FILTER  TIME  CONSTANT  1 

°VEL 

NO.  Of  AtS 

(SEC)  * 

.6400 

.1600 

•713on 

.177oN/At 

1 

.5 

.3392 

.0350 

. 487o„ 

.0469oN/At 

4 

2 

.1800 

.0089 

•3450„ 

.0156oN/At 

10 

5 

.13642 

.0050 

•298°n 

.  0099Ojj/At 

13 

6.5 

.1100 

.0032 

.267cn 

.  0070Cjj/At 

17 

8.5 

.04422 

.0005 

• 167 

.0017oN/At 

44 

22 

Optimum  Filter:  b  ■  a2/(2-a) 


a 

b 

°POS 

°VEL 

FILTER  TIME  CONSTANT 

NO.  Of  AtS 

(SEC)  * 

.6400 

.3010 

.738on 

,342oN/At 

1 

.5 

.3392 

.0690 

.519on 

.093oN/At 

2 

1 

.1800 

.0178 

.3730n 

.031oN/At 

4 

2.0 

.13642 

.0100 

•324on 

.0198oN/At 

5 

in 

• 

CM 

.1100 

.0064 

•29% 

.0l41oN/At 

6 

3 

.04422 

.0010 

•183on 

.0034cN/At 

17 

8.5 

•Assuming  At  ■  0.5  sec 

Oj!  ■  o  of  measured  noise  At  ■  time  interval  between  updates 
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II  FILTER  SIMULATION 

The  digital  simulation  of  the  filters  defined  in  Section  I 
imposed  the  requirement  that  we  also  model  the  ship's  course  as  a 
function  of  time  in  order  to  provide  a  basis  for  computing  the  ideal 
(geometric)  TD  as  well  as  ship's  heading  angle.  It  was  decided  to 
simulate  the  typical  pattern  of  ship  maneuver  in  a  river  channel:  a 
straight-line  segment,  followed  by  a  turn  and  then  another  straight- 
line  segment.  By  repeating  this  pattern  it  would  be  possible,  if 
desired,  to  simulate  the  passage  of  a  ship  over  a  course  of  arbi¬ 
trary  shape  and  length.  Obviously  the  sudden  imposition  and  removal 
of  an  acceleration  normal  to  the  ship's  course  at  the  start  and  at 
the  end  of  the  turn  will  reveal  the  dynamic  response  of  the  filter. 

The  organization  of  the  filter  simulation  program  is  re¬ 
presented  by  the  nine  blocks  in  Fig.  II-l  which  have  some,  but  not 
all,  of  the  properties  of  a  flow  chart.  Each  row  in  the  blocks  of 
the  figure  is  divided  into  two  or  three  columns.  The  number  in  the 
left-hand  column  represents  a  point  of  entry  into  the  block;  the 
number  in  the  right-hand  column  represents  a  point  of  entry  in  the 
same  or  a  different  block  to  which  control  is  passed.  The  center 
column,  when  it  is  present  (blocks  1,  2,  7,  8)  represents  the  set¬ 
ting  of  control  parameters  or  the  testing  of  variables  (generally 
not  explicitly  indicated  on  the  figure)  for  control  purposes. 

Sections  II-l  through  II-9  describe  in  appropriate  detail 
the  operations  performed  in  the  nine  blocks  of  Fig.  II-l.  The 
principal  symbols  used  in  these  sections  are  defined  in  the  follow¬ 
ing  glossary  and  table  of  input  parameters. 
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Glossary  of  Symbols 


x,y 


*,y 


•  •  M 

X/y 

* 

* 

a 

♦ 

T 

TD 

TD 

•• 

TD 

A  A  ^ 

TD,TD,TD 

6TD 

A 

•STD 

TDP 

xp/yp 


Cartesian  position  coordinates  of  the  ship  relative  to 
the  master  LORAN  station. 

First  time  derivatives  of  x,y.  Velocity  coordinates 
of  the  ship. 

Acceleration  coordinates  of  the  ship. 

Ship's  heading  angle,  measured  clockwise  from  north. 

Ship's  turning  rate 

Normal  acceleration  of  the  ship. 

At  any  waypoint,  the  angle  between  the  lines  of  bearing 
to  the  preceding  and  succeeding  waypoints.  (The  ship 
must  turn  through  the  angle  <p  at  each  waypoint.) 


Time  variable. 

Time  difference  between  the  arrival  times  at  the  ship 
of  LORAN  pulses  from  a  slave  and  from  the  master  station 

Time  derivative  of  TD. 

Second  time  derivative  of  TD. 

Filtered  values  of  TD,  TD  an^  TD. 


Error  in  TD.  STD«TD-TD  where  TD  denotes  the  true  value. 

A  A  A 


Error  in  TD.  6TD*TD-TD  where  TD  denotes  the  true  value. 
Value  of  TD  at  an  arbitrary  waypoint,  P 


Cartesian  .position  coordinates  of  an  arbitrary  waypoint, 
P 
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Input  Para*-  iers 


v 

ymax 


dT 


i 


ata 


atk 


atl 


A 

B 

C 

NSIGMA 


Ship's  speed. 

Ship's  maximum  turning  rate. 

Incremental  time  step  in  TO  filter 

Time  delay  in  the  measurement  of  ship's  heading  angle 

Time  delay  in  the  measurement  of  TDs. 

Time  interval  between  TD  measurements  for  TD  filter 

Time  interval  between  measurements  of  ship's  heading  for 
Heading  Filter. 

Gain  in  TD  filter  (position  gain) . 

Gain  in  TD  filter  (rate  gain). 

Gain  in  Heading  Filter 

Standard  deviation  of  noise  in  measured  TDs. 
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II-l  Controller 

The  TD  filter  simulation  consists  of  three  consecutive 
phases:  (1)  a  straight-line  run  of  the  ship;  (2)  a  constant-radius 
turn;  (3)  a  straight-line  run.  Each  phase  requires  the  precompu¬ 
tation  and  setting  of  initial  values  of  certain  parameters.  These 
initializations  are  accomplished  in  Section  II-2,  II-3  and  II-4. 

The  function  of  the  controller  section  is  to  route  con¬ 
trol  of  the  computation  to  the  appropriate  section  as  each  of  the 
phases  is  completed.  At  the  end  of  the  simulation  control  is  passed 
to  Section  II-9  where  the  tapes  for  generating  the  CALCOMP  plots 
of  the  TD  and  TD  errors  are  created. 

11-2  Initialization  (Straight-Line  Segment  Preceding  Turn) 

Input  Parameters:  v,  x^ ,  y^,  Tg 

Output  Parameters:  xq,  yQ,  xq,  yQ,  x,  y, 

TO,  TMAX,  TDi-1,  TD^,  a^,  ^n_1 

At  the  start  of  the  TD  filter  simulation  the  ship  is  on 
a  straight-line,  constant-speed  course  of  duration  Tg  seconds.  This 
is  followed  by  a  turn  at  the  same  constant  speed  and  also  at  a  con¬ 
stant  angular  rate.  The  values  of  i|»,  x  and  y  at  the  start  of  the 
turn  are  precomputed  and  are  denoted  by  x^  and  y^. 

The  initial  conditions  for  the  straight-line  segment  are 

given  by: 

4*0  “  $1  (2.1) 

x_  »  v  siniji  (2.2) 

o  o 
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x 

y 


V  COS^ 


X1  -  X0  TS 


*1  -  TS 
0 
0 


(2.3) 

(2.4) 

(2.5) 

(2.6) 
(2.7) 


The  times  at  the  start  and  at  the  end  of  the  straight- 
line  segment,  TO  and  TMAX  are  given  by 


TO  *  0 
TMAX  =  T, 


(2.8) 

(2.9) 


The  quantities  TDi_1  and  TD^  are  required  in  Section 
II-7  to  initialize  the  TD  filter.  These  are  obtained  from 


TD.  .  =  TD  ,  T  =  -dT. 
i-l  i 


TDi-l  *  TD  ,  T  =  -dTi 


(2.10) 

(2.11) 


in  which  the  TD  and  TD  computations  are  made  in  Section  II-6. 

Initializing  values  required  in  Section  II-8  for  the  Head¬ 
ing  Filter  are  given  by 


an-l  *  0 


Vl  =  *o 


(2.12) 

(2.13) 
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II-3  Initialization  (Constant  Speed  and  Radius  Turn) 


Input  Parameters: 
Input  Variables: 
Output  Parameters: 


v, 


max' 


$ 


x,y,^,T 


xo'  yo'  V  a'  T0'  TMAX 


At  each  way  station  in  the  chain,  the  angle,  <(>,  between 
the  lines  of  bearing  to  the  preceding  and  succeeding  way  stations 
is  the  angle  through  which  the  ship  must  turn  in  the  neighborhood 
of  the  waypoint.  We  assume  the  turn  is  accomplished  symmetrically 
about  the  waypoint  and  at  the  maximum  turning  rate  of  the  ship, 
^max*  Then  f°r  the  turning  segment  the  ship's  turning  rate  is 


^  ^ 


max* 


(3.1) 


The  normal  acceleration  of  the  ship  is 


a 


(3.2) 


The  initial  values  of  x,  y,  and  T  are  obtained  from  current  values 
of  the  simulation  (i'.e.,  at  the  end  of  the  first  straight-line  seg¬ 
ment)  and  these  are  denoted  in 


xo 

*  X 

(3.3) 

yo 

*  y 

(3.4) 

*o 

-  ^ 

(3.5) 

TO 

-  T 

(3.6) 

TMAX,  the  value  of  T  at  the  end  of  the  turn  is  given  by 
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TMAX  a  TO  +  ($>/i  (3.7) 

luaA 


II- 4  Initialization  (Straight-Line  Segment  Following  Turn) 


Input  Parameters: 

v,  Tg 

Input  Variables: 

x,  y,  \p,  T 

Output  Parameters: 

v  yo'  V  V  T0'  TMAX 

At  the  completion  of  the  waypoint  turn,  the  ship  proceeds 
in  a  straight  line  again  for  a  prescribed  period  of  time,  Tg ,  where¬ 
upon  the  simulation  run  is  concluded.  The  initial  conditions  for 
this  segment  are  obtained  from  the  simulation  variables  and  given  by 


x  —  x  (4.1) 

o 

yQ  -  y  (4.2) 

ij>o  =  ip  (4.3) 

TO  =  T  (4.4) 

x„  =  v  sin^  (4.5) 

o  O 

• 

yQ  -  v  cos ipQ  (4.6) 


The  terminating  value  of  time  is  given  by 
TMAX  -  TO  +  Tg 
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II-5  Track  Generator 

Input  Parameters:  xQ,  y  ,  xQ,  yQ,  TO,  ipQ,  i|>,  v,  a 

Input  Variables:  T 

Output  Variables:  x,  y,  x,  y,  x,  y 

The  Track  Generator  computes  the  position,  velocity  and 
acceleration  coordinates  of  the  ship  for  prescribed  values  of  the 
time.  Two  cases  are  considered:  (1)  constant  speed  and  heading 
of  the  ship;  (2)  constant  speed  and  turning  radius  of  the  ship. 
For  the  case  of  constant  speed  and  heading,  the  formulas  for  the 
position  coordinates  are: 


x  “  X-  +  x  (T-TO) 
o  o 

y  -  y0  +  y0(t-to). 


(5.1) 

(5.2) 


In  (5.1)  and  (5.2)  T  denotes  an  arbitrary  time  and  TO  denotes  the 
(initial)  time  at  the  start  of  the  straight-line  run. 

For  a  constant  speed  and  constant  turning  rate  the  formulas 

are: 


♦  +  ty(T-TO) 

(5.3) 

v  sinip 

(5.4) 

v  cosip 

(5.5) 

a  cosip 

(5.6) 

-a  sinip 

(5.7) 

xQ  -  Y/ii>  ♦  (v/ip)  cosipc 

(5.8) 

yc  +  x/i p  -  (v/ip)sinip0 

(5.9) 

EKII-8 
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II-6  TD  Generator 

Input  Parameters;  an,  a12 ,  xp,  yp,  TDP 
Input  Variables:  x,  y,  x,  y 


Output  Varieties: 


TD, 


TD 


The  input  variables  to  the  TD  Generator  are  the  'true1 
position  and  velocity  coordinates  of  the  ship  obtained  from  the 
Track  Generator.  A  TD  (time  difference)  and  its  time  derivative 
are  computed  as  the  linear  combinations  in 

TD  =  TDP  +  a^Cx-Xp)  +  a12(y-yp)  (6.1) 

TD  «  a^  x  +  a12  y.  (6.2) 

The  coefficients  a^  and  a12  in  (6.1)  and  (6.2)  represent 
any  of  the  two  elements  in  the  (two-or  three-rowed)  A-matrix  which 
gives  the  linear  relationship  between  TDs  and  ship's  position  co¬ 
ordinates.  In  other  words,  any  of  the  two  or  three  TDs  may  be  re¬ 
presented  in  the  simulation  by  the  linear  forms  in  (6.1)  and  (6.2) 

The  values'  of  TD  and  TD  given  by  (6.1)  and  (6.2)  represent 
the  geometric  part  of  these  variables  (assuming  noise,  truncation 
errors  and  biases  are  not  present)  for  the  purposes  of  the  simulation. 
They  are  also  used  as^  the  references  or  ’true’  values  for  computing 
the  errors  6TD  and  6TD  defined  in  Section  II-7.  They  are  not  to  be 
confused  with  the  surveyed  TDs  denoted  by  TD^  in  Ref.  131. 


II-7  TD  Filter  with  Dead  Reckoning 


Input  Parameters: 
Input  Variables: 
Output  Variables: 


dTi,  ATa,  ATk,  A,  B,  NSIGMA 


T,  TD,  TD,  TD 

A  *  *  . 

TD,  TD,  STD,  STD 
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Loran  TD  measurements  are  generated  at  intervals  of  ATR 
sec.  A  TD  measurement  available  at  time  is  delayed  or  stale  by 
ATa  sec.  The  TD  measurement  is  derived  from  ten  equally-spaced  TD 
samples,  each  of  which  is  truncated  to  10  nanosec.  The  TD  measure¬ 
ment  is  taken  to  be  the  average  of  the  ten  samples,  also  truncated 
to  10  nanosec.  The  ten  sample  times,  t^,  j-1,..., 10,  are  given  by 

Tj  ■  Tk  -  (ATa+0.45ATr)  +  (j-1)  dT  j  (7.1) 

in  which 

dT .  *  AT  /10 
J  K 

is  the  time  interval  between  TD  samples.  The  time  intervals  ATR, 
ATa#  dTj  and  the  discrete  times  t^,  T^  and  are  illustrated  in 

Fig.  II-2 . 

Let  TDj  denote  a  true  value  of  a  TD  at  time  and  XRAN^ 
the  error  in  the  TD  sample.  The  TD  sample,  truncated  to  10  nanosec 
(with  the  5.0  nanosec  truncation  bias  removed),  is  denoted  by  TDS ^ 
and  given  by 

TDS j  «  10 [0. 1 (TD j+XRAN j ) ]  +  5.0  (7.2) 

in  which  lx]  denotes  the  integer  part  of  x. 

Finally,  the  TD  measurement,  which  is  the  mean  of  the  ten 
samples,  truncated  to  10  nanosec,  is  denoted  by  TDM  and  given  by 

10 

TDM.  »  10(0.1(1/10)  I  TDS.]  +  5.0.  (7.3) 

K  j-1  3 

In  the  filter  simulation  the  TD  error  XRAN  was  represented  by  a 
normally  distributed  random  variable  with  zero  mean  and  standard 
deviation  denoted  by  NSIGMA. 
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The  integration  time  step  in  the  TO  filter  is  denoted  by 
dT^.  During  intervals  between  TD  measurements  dead-reckoning  equa¬ 
tions  arc  used  to  produce  the  filter  outputs,  denoted  by  TD^  and 
TD^,  every  d^  seconds.  The  dead- reckoning  equations  are: 


TD^  ■ 

TDi-l  + 

(TDi -i+0. 5  TD  dTiidTi 

(7.4) 

A 

A 

A 

• 

• 

•• 

(7.5) 

TDi  * 

TDi-l  + 

TD  dTi 

*1-1 

II 

o> 

(7.6) 

A 

=  TD. 
l 

(7.7) 

/\ 

•• 

The  value  of  TD  used  in  (7.4)  and  (7.5)  is  taken  from  the  current 
output  of  the  Heading  Filter. 

When  a  new  TD  measurement  first  becomes  available  a  differ- 

A  A 

ent  set  of  filter  equations  are  used  to  modify  TD^  and  TD^.  In  this 
circumstance  the  seven  applicable  filter  equations  are: 


TD^  »  TDi 

A  A 

TD^  »  TDi 

A  A  A 

Zk  -  TDM^  +  ATA(TDk+0.5  AT&  TD)  -  TD^ 

A  A 

TDk  »  TDk  +  (A+BUTA/ATR)  )Zk 

A  A 

TDk  -  TDk  +  (B/ATR)Zk 
TD^  ■  TDk 

A  A 

TDi  -  TDk 


(7.8) 

(7.9) 

(7.10) 

(7.11) 

(7.12) 

(7.13) 

(7.14) 
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Finally,  the  errors  in  the  outputs  of  the  TD  Filter  are 


given  by 


STD  -  TDa  -  TD 


(7.15) 


STD  ■  TDa  -  TD 


(7.16) 


II- 8  Heading  Filter 

Input  Parameters:  ATL,  C,  a^,  a12,  v 

Input  Variables:  ip,  x,  y 

^  A 

Output  Variables:  \p  ,  TD»  STD 

The  Heading  Filter  uses  a  measurement  of  ship's  heading 
angle  to  produce  a  filtered  value  of  ship's  normal  acceleration 

A 

from  which  is  obtained  TD  used  in  the  dead-reckoning  computations 
in  (7.4),  (7.5)  and  (7.10)  of  the  TD  Filter.  The  time  interval 
between  heading  measurements  (and  hence  updates  in  the  value  of 
TD)  is  denoted  by  ATL  and  the  measurement  of  ship's  heading  angle 
is  delayed  or  stale  by  dT^. 

The  true  value  of  ship's  heading  angle  at  time  t=t  is 

n 

denoted  by  ipn  and  expressed  in  radians.  To  express  the  heading  angle 
in  degrees  we  write 


r|>n(180/w)  . 


(8.1) 


The  measured  \p  is  expressed  in  degrees  and  truncated  to  the  nearest 

whole  degree.  We  convert  this  measured  value  back  to  radians  and 

denote  its  value  by  tbt.  We  have 

n 
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% 


li|£)  (tt/180) 


(8.2) 


and 


^n-l  “  ^n-l1  <1T/180) 


(8.3) 


where  ^  .  denotes  the  previously  measured  value  of  ij;  at  t*t  -  AT_  . 

n  L 

Let 


44,t  *  *n  *  *n-l' 


(8.4) 


The  filtered  value  of  ship's  normal  acceleration,  denoted  by  a,  is 
then  given  by 


*n  "  Sn-1  + 


(8.5) 


The  components  of  an  are 


x  »  a  cos  ill 
n  n  n 


yn  -  -an  sin  V 


(8.6) 

(8.7) 


Then  the  filtered  value  of  TD  is  obtained  from 


TD  *  *n  all  +  a12 


(8.8) 


and 


TD  ■  x  a^  +  y  a^2* 

A 

The  error  in  TD  is  given  by 


(8.9) 
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A 

5TD  -  TD  -  TD. 


(8.10) 


Finally, 


i  ®  a  (8.11) 

n-i  n 

Vi  -  v  (8-12> 

II-9  CALCOMP  Plots 

At  the  conclusion  of  the  simulation  control  is  passed  to 
block  9  in  which  the  tapes  for  creating  the  CALCOMP  plots  are  made. 
In  this  block  the  means  and  standard  deviations  of  the  TD  error  and 
its  derivative  are  also  computed. 

After  a  plot  tape  is  made  the  program  reads  the  input  tape 
to  see  if  there  is  another  simulation  to  be  made. 
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U.S.C.G  LOR AN -C  TD  FILTER  SIMULATION 


START 


2 


INITIALIZE  (1) 

50 

LGO-4 

100 

90 

300 

I 


CONTROLLER 

14 

KGOsO 

30 

KGO*  1 

50 

30 

KG0*2 

60 

KG0S3 

70 

KG0«4 

19 

3 


INITIALIZE  (2) 

60 

80 

4 


INITIALIZE  (3) 

70 

80 

80 

300 

5 


TRACK  GENERATOR 

100 

200 

7 


T 

D  FILTER 

300 

30 

400 

315 

LGO-I 

100 

305 

LG0*3 

316 

300 

6 


TD  GENERATOR 

LGO*i 

315 

200 

LG0«2 

420 

LG0«3 

316 

LG0«4 

90 

9 

CALCOMP  PLOTS 

l|9_ 

END 

8 


HEADING  FILTER 

400 

305 

410 

410 

LG0=2 

100 

420 

305 
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III  FILTER  PERFORMANCE 

The  waypoint  chosen  for  the  sample  simulations  presented 
in  this  section  was  numbered  2  in  the  numbering  scheme  employed  in 
the  spring  of  1979.  This  was  the  second  waypoint  on  the  downbound 
summer  route  and  required  a  turn  through  64.5  deg. 

Numerical  values  of  system  parameters  common  to  all  simu¬ 
lations  reported  here  are  given  in  Table  III-l.  The  turn,  executed 
at  0.5  deg/sec,  required  129.0.  sec.  Straight-line  runs  of  10  sec 
preceded  and  followed  the  turn.  Therefore  a  simulation  run  extended 
over  an  interval  of  149  sec. 

In  Plots  1,  2,  3  and  4  the  gain  C  in  the  Heading  Filter 
is  given  the  values  0,  0.2,  0.5,  0.9,  respectively.  The  means  and 
standard  deviations  of  the  TD  error  and  its  derivative  for  these 
four  simulations  are  tabulated  in  Table  III-2. 

In  each  of  the  four  plots  the  scale  for  the  TD  error  in 
nanosec  is  given  at  the  left  of  the  plot  and  the  TD  error  scale  in 
nanosec/sec  is  on  the  right.  Observe  the  changes  in  scale  which 
occur  in  going  from  plot  to  plot. 

In  each  plot  the  lower  curve,  which  also  has  the  larger- 
amplitude,  high-frequency  oscillations,  represents  the  TD  error. 

In  Plot  #1  the  gain  C  in  the  heading  filter  is  zero,  so 
that  this  run  corresponds  to  the  conventional  second-order  filter 
without  acceleration-aiding.  Note  the  large  velocity  lag  (TDDT 
error)  which  exists  at  the  end  of  the  run.  The  TD  error  is  nearly 
240  nanosec  -  a  value  much  too  large  to  be  tolerated. 

In  Plots  #2,  3,  and  4  as  the  gain  C  is  increased  the  TD 
error  is  correspondingly  decreased. 

While  the  runs  shown  in  the  plots  represent  the  case  of 
a  critically  damped  filter  (compare  values  of  A  and  B  in  Table  III-l 
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with  a  and  b  in  Table  1-1)  we  do  not  intend  to  imply  that  this  is 
the  'best1 choice  of  filter  9ains  for  an  operational  system.  They 
are  presented  because  they  indeed  exhibit  the  behavior  one  might 
expect  of  a  critically  damped  filter  and  thus  add  to  our  confidence 

in  the  simulation. 
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Table  III-l 

Numerical  Values  of  System  Parameters 


V 

= 

10.0  knots 

^max 

a 

0.5  deg/sec 

dTi 

a 

0.01  sec 

dT^ 

= 

0.5  sec 

ata 

a 

0.5  sec 

ATk 

a 

1.0  sec 

AT 

L 

= 

1.0  sec 

A 

a 

0.04422 

B 

a 

0.0005 

NSIGMA 

a 

5.0  nanosec 

d-iii-3 


E 
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TO:  C.  R.  Edwards 

FROM:  V.  Schwab 

SUBJECT:  The  Effect  of  Uncompensated  TD  Truncation  Biases 

on  TD  and  T*D  Errors  in  the  USCG  LORAN-C  HHE 
Navigation  Program 

REFERENCE:  S2R-79-166,  Design,  Simulation  and  Performance  of 

Digital  TD  and  Heading  Filter  for  the  USCG  LORAN-C 
HHE  Navigation  Program,  by  L.  J.  Levy  and  V.  Schwab, 
October  22,  1979 

In  the  USCG  LORAN-C  HHE  Navigation  Program  a  10-nanosec 
truncation  is  performed  on  each  of  ten  equally-spaced  TD  samples 
which  are  then  averaged  in  the  receiver.  The  average  of  the  ten 
truncated  samples  is  itself  then  truncated  at  the  10-nanosec  level 
to  produce  a  TD  measurement.  The  truncation  of  the  samples  pro¬ 
duces  a  5.0-nanosec  bias  in  the  average  and  the  truncation  of  the 
average  creates  an  additional  5.0  nanosec  bias  in  the  TD  measure¬ 
ment  for  a  total  bias  of  10  nanosec. 

It  is  desirable  to  compensate  for  these  biases  and  quite 
easy  to  do  so.  A  compensation  of  5.0  nanosec  should  be  added 
either  to  each  TD  sample  in  the  receiver  or  to  the  average  of  the 
ten  TD  samples  before  it  is  truncated.  Then  the  truncated  average 
should  be  compensated  by  the  addition  of  5.0  nanosec.  However  if 
separate  5.0-nanosec  compensations  are  not  made  to  the  samples  and 
the  measurement,  a  single  10.0-nanosec  compensation  may  be  added  to 
the  measurement  with  almost  equal  beneficial  effect. 
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IBM  3033  Computer  Simulations  Which  Show  the  Effect  of  Compensated 
and  Uncompensated  Truncation  Biases 

The  DSC6  TD  Filter  Simulation  Program  described  in  the 
reference  was  used  to  study  the  behavior  of  TD  and  TD  errors  in 
the  presence  of  biases,  uncompensated  and  compensated,  resulting 
from  TD  sample  and  measurement  truncations. 

In  Table  1  the  results  of  five  simulation  runs  are  shown. 
The  filter  gains  and  system  constants,  defined  in  the  reference 
memorandum,  had  the  same  values  indicated  in  the  Table  in  each  of 
the  five  runs.  What  was  varied  from  run  to  run  were  the  compensa¬ 
tions  made  for  the  biases  resulting  from  the  sample  and  measurement 
truncations. 

In  Run  1  no  compensation  is  made  for  truncation  bias  and 
the  mean  of  the  TD  error  is  11.268  nanosec  -  close  to  the  total 
10.0-nanosec  bias  present  in  the  TD  measurement. 

In  Runs  2  and  3  the  sample  and  measurement  truncations 
are  compensated,  but  not  both  in  the  same  run.  As  a  result  the 
mean  TD  errors  are  not  far  from  the  expected  5.0  nanosec  bias. 

The  proper  compensations  for  the  truncation  biases  are 
made  in  Run  4  and  here  the  mean  TD  error  is  less  than  a  quarter  of 
a  nanosec. 

The  single  10.0-nanosec  compensation  made  in  Run  5  is 
not  quite  as  effective  as  the  dual  compensations  of  Run  4;  still, 
the  mean  TD  error  in  this  case  is  only  a  little  greater  than  one 
nanosec. 

The  simulations  in  Table  1  were  made  for  a  stationary  ship 
(V»0)  and  no  heading  filter  was  used  (C«0) .  The  filter  was  slightly 
underdamped  (B>0.005). 
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TO: 

FROM: 

SUBJECT: 

REFERENCES: 


C.  R.  Edwards 

L.  J.  Levy/V.  Schwab 

Formulas  for  Making  Short-Range  Predictions  of 
Ship's  Course  in  the  USCG  LORAN-C  HHE  Navigation 
Program 

[1]  APL/JHU  S2R-79-162,  Methods  for  Obtaining 
Approximate  and  Exact  Solutions  to  a  LORAN 
Navigation  Equation,  by  L.  J.  Levy  and 
V.  Schwab,  October  8,  1979 

12]  APL/JHU  S2R-79-166,  Design,  Simulation  and 
Performance  of  Digital  TD  and  Heading 
Filters  for  the  USCG  LORAN-C  HHE  Navigation 
Program,  by  L.  J.  Levy  and  V.  Schwab,  October 
22,  1979 


It  may  prove  helpful  to  a  pilot  in  steering  his  ship  to 
provide  him  with  pictures  (CRT  displays)  of  the  predicted  position 
of  the  ship  a  short  time  into  the  future.  By  observing  how  the 
predicted  position  responds  to  changes  in  engine  speed  and  rudder 
position  he  may  be  significantly  aided  in  anticipating  the  maneuvers 
required  to  follow  a  prescribed  or  desired  course. 

The  purpose  of  this  memorandum  is  to  provide  the • formulas 
for  making  simple  short-range  predictions  of  ship  position  based  on: 
present  ship's  velocity  alone;  present  ship's  velocity  as  well  as 
present  turning  rate.  We  consider  how  the  initial  conditions  of  posi¬ 
tion  and  velocity  might  be  derived  from  LORAN-C  measurements.  It  is 
proposed  that  initial  turning  rate  be  obtained  from  filtered  dif¬ 
ferences  of  gyrocompass  readings. 


1.  Formulas  for  Predicting  Ship's  Course 
We  define  the  following  terms: 
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x0/Yq  present  cartesian  position  coordinates  of  the  ship 

Xq'^O  present  velocity  coordinates  of  the  ship 

present  heading  angle  of  the  ship  (measured  clock- 
wise  from  north) 

• 

\|i0  present  turning  rate  of  the  ship 

vQ  present  speed  of  the  ship 


T  time  interval  over  which  course  predictions  are  to 
be  made 


D 

At 

t(k) 

Ax(k) ,Ay (k) 


distance  interval  over  which  course  predictions  are 
to  be  made 

time  interval  between  predictions 

time  of  kth  prediction  (k=*0 , 1, 2  , . . .  ,KMAX) 

predicted  coordinates  of  the  ship  at  time  t(k)  rela 
tive  to  predicted  position  at  time  t(k-l) . 


Initially  (t (0) =0)  the  ship  is  at  x=xQ  and  y=yQ.  At  t(l)  = 
At  the  predicted  coordinates  are  x1=xQ+Ax(l)  and  y^yg+Ayd)  ;  at 
t(2)=2At,  X2=xq+Ax(1)  +  Ax(2)  and  y2=  y^+Ayd)  +  Ay(2)  .  We  shall 
compute  the  increments  Ax(k)  and  Ay(k)  rather  than  x^  and  y^  since 
the  CRT  plotter  makes  use  of  these  increments.  We  consider  two  cases 
one  in  which  the  prediction  is  based  upon  the  ship's  present  trans¬ 
lational  velocity  alone  and  the  other  in  which  account  is  also  taken 
of  the  present  turning  rate  of  the  ship. 

Given:  xQ ,yQ ,yQ ,T  or  D 


Compute; 


At=T/KMAX  or  At=D/ (KMAX  VQ) 
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DO  1  K=1 ,  KMAX 
t(K)  =  KAt 
Ax(K)  =  xQAt 
1  Ay (K)  =  yQAt 

X0 ' ^0 ' xo ' ^0 ' ^0 ,  ^  or  ® 

At  =  T/KMAX  or  At  =  D/ (KMAX  VQ) 
A«P  =  ipQAt 


sin  *0  -  *0/vo 
cos  =  y0/vo 

sin  Aip  =  Aip  -  Aip3/3 1  +  ... 
cos  Aip  =  1  -  Aip2/2  J  +  ... 

DO  2  K=1 ,  KMAX 

sin  ip  =  sin  ip^  cos  Aip  +  cos  ip^  sin  Aip 
cos  ip  =  cos  ipQ  cos  Aip  -  sin  ipg  sin  Aip 
t (K)  *  KAt 

Ax(K)  =  (Vq/iPq)  [cos  ipQ  -  cos  ip] 
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Ay(K)  =  (v0/ij)0)  [sin  -  sin  <Ji] 
sin  =  sin  ^ 

2  cos  =  cos  ip 

2.  Source  of  xQ  and  yQ 

The  present  horizontal  coordinates  of  the  ship,  xQ  and 
y^,  are  obtained  as  the  exact  solution 


of  equation  (4.1)  of  Ref.  [1],  The  iterative  scheme  for  obtaining 
Zq  is  given  by  (4.12). 

3.  Source  of  xQ  and  yQ 

It  would  be  possible  to  produce  x.  and  y  as  the  exact 

•  v  U 

solution  Zq,  of  a  LORAN  rate  equation  analogous  to  (4.]).  For  the 
intended  purpose  that  accuracy  is  not  required  and  the  computational 
burden  could  not  be  justified.  On  the  other  hand  the  linear  G-matrix 
approach  would  probably  result  in  errors  that  are  larger  than  we  can 
stand. 

Instead  we  recommend  that  the  x-  and  y-  outputs  of  the 
LORAN  equation  be  differentiated  and  then  passed  through  a  first- 
order  filter  to  produce  x0  and  y0>  In  Ref.  [2]  the  time  interval 
between  TD  updates  is  denoted  by  AT„.  We  shall  assume  here  that 
predictions  based  on  the  updated  values  of  x^  and  y^  are  to  be  made 
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at  the  same  time  frequency.  We  denote  the  differentiated  LORAN 
position  outputs,  xQ  and  yQ,  by  x*  and  y*  and  write 

.*  xolt>  '  xo(t  '  atk> 

x°  =  Atk 

...  V0(t)  -  y0(t  " 

yo  StJJ 

The  filtered  velocity  coordinates  xQ(t)  and  yQ (t)  are  then  obtained 
from  the  first-order  filter  equations 


xQ(t)  =  xQ(t  -  ATr)  +  H[xJ  -  xQ(t  -  ATRn 

yQ(t)  =  yQ(t  -  ATr)  +  H{yJ  -  yQ(t  -  ATR)  ] 

in  which  the  gain  constant  H  will  be  chosen  to  produce  the  desired 
performance  of  the  predictor.  H  will  depend  on  the  gain  settings 
in  the  TD  filter  which  control  the  noise  levels  in  xQ  and  y^. 


4.  Source  of  i 
_ o 

The  first-order  heading  filter  which  supplies  a  filtered 
value  of  TD  is  defined  by  equation  (8.5)  of  Section  II  of  Ref.  [2] 
and  is  reproduced  below: 

an  *  Sn-1  +  C[v(A*VATL)  -  an_^. 
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In  the  equation  above  an  =  is  the  filtered  acceleration  normal 
to  the  path  of  the  ship,  C  is  the  gain  constant,  v  denotes  ship's 

f  »  •  A  A  f 

speed,  AT_  is  the  time  interval  between  a  and  a  ,  and  Atf/  denotes 
L  n  n-1  r 

the  difference  of  truncated  gyrocompass  readings. 

One  possibility  for  obtaining  ^  for  the  predictor  would 
be  to  use  the  formula 


♦o  - 


Vv0 


The  possible  disadvantage  resulting  from  this  source  for  iji  is 
the  value  of  C  used  in  the  a-filter.  For  the  TD-filter  a  large 
value  of  C  (near  unity)  is  to  be  preferred  in  order  to  obtain  a 
rapid  response  to  a  turn.  For  the  predictor  application,  a  smaller 
value  of  C  might  prove  to  be  a  better  choice.  Therefore  we  recom¬ 
mend  that  the  following  separate  ^-filter  be  employed  for  the  pre¬ 
dictor: 

5>0(t)  =  ^Q(t-ATL)  +  K(A^VATL  -  $>0(t-ATL)]. 

In  the  equation  above  K  denotes  the  gain  constant  and  A^t  is  the 
same  term  which  appears  in  the  heading  filter. 
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November  26,  1979 


TO: 


C.  R.  Edwards 


FROM:  V.  Schwab 

SUBJECT:  Additional  Performance  Evaluations  for  the  Second- 

Order  a,b-Filter 

REFERENCE:  S2R-79-166  Design,  Simulation  and  Performance  of 

Digital  TD  and  Heading  Filters  for  the  USCG  LORAN-C 
HHE  Navigation  Program,  by  L.  J.  Levy  and  V.  Schwab, 
October  22,  1979 


In  Table  1-1  of  the  reference,  reproduced  in  this  memo¬ 
randum,  we  present  evaluations  of  some  of  the  performance  charac¬ 
teristics  of  the  a,b-filter  used  to  filter  LORAN  TDs  in  the  USCG 
HHE  Navigation  Program.  Table  1-1  offers  examples  of  two  filter 

designs:  (1)  a  critically-damped  filter  defined  by  the  gain  rela- 
2 

tion,  4b=(a+b)  ;  (2)  an  optimum  filter,  defined  by  the  gain  rela- 

2 

tion,  b=a  /(2-a) ,  which  results  in  a  slightly  underdamped  filter. 

We  have  added  to  Table  1-1  a  third  case  in  which  the  rate 
gain  b  is  chosen  to  be  the  mean  of  the  values  of  the  first  two 
cases,  i.e.,  b=0.5{2-a  -  2/1-a  +  a  /(2-a)}. 
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Table  1-1 

Summary  of  Some  Characteristics  of 
The  Second-Order  Filter 


Critically  Damped  Filter:  4b 


(a+b)2 


a 

b 

°POS 

FILTER  TIME 

CONSTANT 

°VEL 

No.  of  Ats 

(SEC)  * 

.6400 

.1600 

*  713°n 

•177oN/At 

1 

.5 

.3392 

.0350 

.487on 

.0469oN/At 

4 

2 

.1800 

.0089 

.345on 

.0156oN/At 

10 

5 

.13642 

.0050 

.298on 

.0099aN/At 

13 

6.5 

.1100 

.0032 

.  2  6  7  o 

N 

. 0070o  /At 

17 

8.5 

.04422 

.0005 

.!67cn 

. 0017cN/At 

44 

22 

Optimum  Filter:  b  =  a2/(2-a) 


a 

°VEL 

FILTER  TIME  CONSTANT 

No.  of  Ats 

(SEC)  * 

.6400 

.3010 

.738oN 

.  342a_,/At 

•5 

.3392 

.0690 

’  33"^°N 

. 093a„/At 

2 

1 

.1800 

.0178 

*373on 

. 0  31o../At 

4 

2.0 

.13642 

.0100 

*324on 

.0198oN/At 

5 

2.5 

.1100 

.0064 

.290on 

. 0l41aN/At 

6 

3 

.04422 

.0010 

•  183o 

N 

.0034aN/At 

17 

8.5 

♦Assuming  At  =  0.5  sec 

a.,  =  a  of  measured  noise 
N 
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At  =  time  interval  between  updates 
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December  18,  1979 


TO: 


C.  R.  Edwards 


FROM:  L.  J.  Levy/V.  Schwab 

SUBJECT:  Exact  Method  for  Obtaining  Horizontal  components  of 

Ship’s  Velocity  from  TD  Measurements 

REFERENCES:  [1]  S2R-79-187,  Formulas  for  Making  Short-Range  Pre¬ 

dictions  of  Ship's  Course  in  the  USCG  LORAN-C 
HHE  Navigation  Program,  by  L.  J.  Levy  and 
V.  Schwab,  November  20,  1979 

[2]  S2R-79-162,  Methods  for  Obtaining  Approximate 
and  Exact  Solutions  to  a  LORAN  Navigation  Equa¬ 
tion,  by  L.  J.  Levy  and  V.  Schwab,  October  8, 

TSTT 

[3]  S2R-79-008,  G-Matrix  Analysis  for  USCG  LORAN -C 
HHE  Navigation  Program,  by  L.  J.  Levy  and 

V.  Schwab,  January  19,  1979 


In  Ref.  [1]  we  proposed  that  the  horizontal  velocity  com¬ 
ponents  of  the  ship  be  determined  by  passing  differentiated  navi¬ 
gated  positions  through  first-order  filters.  The  alternative  method, 
much  to  be  preferred  if  the  computing  capacity  is  available,  is  to 
make  use  of  the  relation 


5q- 


gq 


in  which  TD-  is  the  vector  of  TD  time  derivatives  obtained  from  the 
— Q 

LORAN  TD  filter  and  Gq  is  the  G-matrix  evaluated  at  Z=Zq,  the  current 
navigated  position  of  the  ship.  In  this  method,  the  velocity  compu¬ 
tation  is  performed  after  the  position  fix  has  been  obtained  (i.e., 
after  the  value  of  has  been  determined  by  means  of  the  iterative 

procedure  outlined  in  Section  4  of  Ref.  [2]).  The  relation  given 
•  • 

above  between  Z_  and  TD„  is  exact  since  Gn  is  evaluated  at  Z=Z_.  It 
— y  — -q  u  —  ~v 
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therefore  fully  exploits  the  velocity  information  contained  in  the 
TPs  that  are  generated  by  the  TD  filter. 

In  the  following  sections  we  list  the  computations  required 
to  produce  Z^  for  the  cases  where  there  are  three  and  two  TD  measure¬ 
ments  available. 

1*  Velocity  Computation  Employing  Three  TPs 

In  the  notation  of  Ref.  [3],  the  linearized  relation  be¬ 
tween  horizontal  position,  represented  by  the  2x1  vector  and 
LORAN  time  differences,  represented  by  the  3x1  vector  TD,  is  given  by 

TD_  -  TD  *  A  (Z_  -  Z  )  .  (1.1) 

— Q  — p  P  —Q  — p 

In  (1.1)  the  elements  of  the  3x2  A  -matrix  depend  upon  Z^  and  the 

p  ~ P 

coordinates  of  the  master  and  slave  stations.  If  we  'solve'  (1.1) 

for  (Z_-Z  )  we  obtain 
— Q  — P 


*0 


-  z  =  G  (TD_  -  TD  ) 
— p  p  — -Q  — P 


(1.2) 


in  which  G= (A  A* ) _1A*- 
P  P  P  P 

Equations  (1.1)  and  (1.2)  are  linear  approximations  in 
which  p  denotes  a  fixed  waypoint  and  Q  a  variable  nearby  point.  If 
we  differentiate  (1.2)  with  respect  to  the  time  we  obtain  the  ap¬ 
proximation 


GP^Q* 


However,  at  p*»Q  equation  (1.3)  becomes  exact,  i.e., 


(1.3) 


(1.4) 
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and  provides  the  basis  for  the  computation  of  ship's  velocity.  In 
Section  3  we  present  a  comparison  of  numerical  results  obtained  by 
means  of  (1.3)  and  (1.4). 

We  next  list  the  equations  used  to  compute  the  elements  of 
the  Gq  matrix.  First  let 


h  - 


(1.5) 


represent  the  navigated  horizontal  position  coordinates  of  the  ship. 
Then  r ^  denotes  the  slant  range  from  the  point  Q  to  the  master  sta- 

'm'-‘m'  “2j 


tion  (x  ,y  )  and  r, .  the  slant  range  from  the  point  Q  to  the  jth 


slave  station  (x^,y^).  Evidently 

r»  ■  ["q  - 


*-[ 


2  2~1  ^ 

V  *  %- V  J 


1/2 


(Xf*  •  x.)  +  (yn  -  y.)  L  jslf2f3* 


] 


The  elements  of  the  A^-matrix  are 


AQjl  -  ij 

|  -Q 

L 

AQj2  -  i 

1^ 

V 

L 

ih  -  -  xq  -  V)  jMl2. 

vL  r2j  r2  J  -  3  ^ 


-  y-s  -  yn  - 


,  j  -  1»2,3 


where  v  denotes  the  signal  propagation  velocity. 
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Next  let 


Dll  ■  AQll2  +  AQ212  +  AQ312 

D12  -  AQ11*AQ12  +  AQ21*AQ22  +  AQ31*AQ32 

D21  »  D12 

D22  =  AQ122  +  AQ222  +  AQ322 
A  a  D22*Dll  -  D21*D12. 

Then 

GQlj  -  (D22*AQ jl  -  D21*AQj2)/A,  j«l,2,3 

GQ2 j  -  (Dll*AQj2  -  Dl2*AQjl)/A,  j«l,2,3. 

Finally 

XQ  -  GQ11*TD1  +  GQ12*TD2  +  GQ13  TD3 
yQ  -  GQ21*TD1  +  GQ22*TD2  +  GQ13  TD3. 

2.  Velocity  Computation  Employing  Two  TDe 

e  e  • 

In  Section  1  It  was  assumed  that  three  TDs,  TDl ,  TD2  and 

•  e  • 

TD3  were  available  to  determine  and  y^.  Here  we  shall  assume 
that  only  two  are  available  and  these  two  we  shall  denote  by  TDj 
and  TDk.  The  possible  values  for  j  and  k  are  the  pairs: 

j  •  1  ,  *  -  2 

j  -  1  ,  *  -  3 

j  -  2  ,  k  -  3 


No  other  values  of  j  and  k  are  admitted 
AQ  Is  the  2x2  matrix 


f 
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AQ 


~AQjl 

AQkl 


AQj2~ 

AQk2 


whose  elements  are  defined  by  expressions  given  in  Section  1.  Let 
A  »  AQjl*AQk2  -  AQkl*AQj2 

and 

GQll  -  AQk2/A 
GQ21  -  -AQkl/A 
GQ12  -  -AQj2/A 
GQ22  -  AQj2/A. 

Then 

Xq  -  GQll*TDj  +  GQl2*TDk 

•  •  • 
yQ  »  GQ21*TDj  +  GQ22*TDk 


3.  Comparison  of  Numerical  Results  Obtained  Using  G  and  G  Matrix 

-  —  p  SL 


In  this  section  we  compare  numerical  results  computed  for 

the  St.  Mary's  River  LORAN-C  mini  chain  using  equations  (1.3)  and 

(1.4)  of  Section  1.  The  approximation  given  by  (1.3)  which  uses  the 

G  matrix  would  offer  some  advantages  in  that  the  G_  matrices  are 
P  P 

precomputed  and  stored  on  tape.  However,  as  we  shall  see,  the  errors 

in  the  velocity  determined  by  means  of  the  G^  matrix  can  be  quite 

large. 


H-5 


THf  JOHNS  HOPKINS  UNI  VCRS)  TV 

APPLIED  PHYSICS  LABORATORY 
laurel  Maryland 


S2R-79-200 
Page  -6- 


In  Table  I  we  list  the  along-track  (ATE)  and  cross-track 
(CTE)  navigation  errors  in  position  and  velocity  at  a  number  of 
points  in  the  St.  Mary's  River.  The  index  I  denotes  the  waypoint 
number  and  K»1  means  downbound,  summer.  S  denotes  linear  distance 
between  successive  waypoints  and  ITER  represents  the  number  of  itera 
tions  of  the  navigation  computation.  The  navigation  position  errors 
are  expressed  in  feet. 

The  navigation  velocity  errors  in  Table  I  are  zero  since 
the  Gg-matrix  is  used  to  compute  Z_.  The  components  of  ship's  velo¬ 
city  remain  constant  between  successive  waypoints  and  ship's  speed 
is  12  mph. 

Table  II  is  the  same  as  Table  I  except  that  the  G  -matrix 
•  P 

is  used  to  compute  Z.  In  this  case  the  along-track  and  cross-track 

velocity  errors  are  no  longer  zero  but  increase  monotonically  with  S 
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E*M.T  SOLUTIONS  CF  LUKAi*  cWOAT lLNS  FUi?  HCAl2CNT<t  POSITION  ANu  VELOCITY  COOPCINATES 


l  k 

1  1 

1  1 

l  1 

l  1 

l  L 

1  1 

l  l 

A  i 

1  1 

1  I 

1  1 

J  1 

2  k 

1  I 

2  i. 

2  A 

A 
A 
A 

I 
A 

A 
A 


A 
i 
A 
A 
A 
A 

6  A 
o  I 
«  1 
6  A 

6  k 

7  A 

V  A 
7  1 

7  A 
7  A 

6  A 

•  A 

9  A 
9  A 
9  A 
9  A 


i  (kYOi) 

If.k 

xTc  1  FT ) 

CTi  |P1> 

1  It* 

ATP  (6®r1) 

lip  inpnl 

At  01  IMP*) 

YD01  (“PH  I 

0.0 

A 

0.0 

0.0 

n 

O.OCCOO 

c.ooooo 

7.95 

— F  .46 

A. AO 

5 

0.05 

-0.01 

0 

O.OOOoO 

0. coooo 

7.V6 

-6.46 

2.20 

5 

0.03 

-O.tO 

0 

c.occco 

-C.COCCU 

7.95 

-6.58 

A.  j  0 

6 

-0.1 1 

0.01 

0 

O.OOO'iO 

-C.OOOGO 

7.95 

-».V4 

4.  AO 

7 

0.26 

-0.03 

0 

O.OC. '<00 

-O.COOOO 

7.V» 

-8.5P 

9.3c 

'» 

0.2. 

-O.C2 

0 

0.00000 

0.00000 

7.91 

—8  .56 

6.60 

AA 

0.2LI 

-0,o2 

0 

0.00000 

-c.ooooo 

7.95 

-6.58 

7.70 

15 

0.25 

-u.02 

c 

e.oooco 

0.000 00 

7.95 

—6.  “8 

e.eo 

15 

0.20 

-0.03 

0 

-O.Ul’CCO 

-o.eoooc 

7.95 

-8.96 

4.90 

A  7 

0.3  / 

-0.05 

0 

c.ooooo 

-C.OOOOO 

7.95 

-4.5* 

U.oo 

19 

0.*-9 

*0.0*1 

0 

-0.00)09 

-0. coo 00 

7.95 

—8. 98 

11.14 

<0 

-0.37 

0.03 

0 

C.OCC'k.O 

-o.coooo 

7.91 

-6.96 

0.0 

A 

0.0 

0.0 

0 

o.ooonu 

-0.00000 

11.93 

3.22 

l.GO 

3 

O.OC 

o.co 

c 

c.ooooo 

-o.coooo 

11.53 

3.22 

2.00 

-* 

-0.01 

-O  .01 

0 

O.OOuOO 

0.00000 

11.53 

3.32 

3.00 

5 

0.02 

0.01 

0 

0.00000 

-0.00000 

11.53 

3.32 

A. 00 

6 

-0.03 

-0.C2 

0 

-0.0000  0 

-o.coooo 

11.13 

3.22 

9.00 

7 

0.04 

0.03 

0 

-C.OCGOO 

o.coooo 

11.53 

2.32 

5.42 

7 

0.08 

0.05 

0 

-O.OOOCO 

-0.00000 

11.32 

3.32 

0.0 

1 

-0.00 

-0.00 

0 

-0.0000  0 

o.coooo 

10.64 

5.56 

0.56 

2 

-0.02 

-0.02 

0 

-o.oooeo 

O.bOOOO 

10.64 

5.56 

0.0 

A 

-0.00 

-O.CO 

0 

—O.OCCOO 

-O.OOOCO 

5.27 

10.78 

A. 00 

2 

0.00 

0.00 

0 

-0.00000 

-0.00000 

5.27 

10.76 

1.56 

2 

0.01 

0.02 

0 

-O.OOOUO 

-0.00000 

5.27 

10.76 

0.0 

A 

-O.OC 

-0.00 

0 

-0.00000 

-O.COOOO 

9.63 

7.18 

1.00 

3 

O.OC 

o.co 

0 

-O.OOO0O 

-O.OOOCO 

9.63 

7.16 

>2.00 

3 

0.03 

0.04 

0 

-O.OOOCO 

-o.coooo 

9.63 

7.14 

3.00 

4 

-0.02 

-0.03 

0 

-O.OOOUO 

-e.oooco 

4.63 

7.16 

4.00 

5 

0.02 

0.02 

0 

-o.ooot-o 

-o.coooo 

*.63 

7.16 

4.29 

5 

0.03 

0.04 

0 

C.OOCO  0 

c.coooc 

9.63 

7.16 

0.0 

1 

-o.oo 

-0.00 

0 

0.00000 

0.00000 

1}  .67 

2.79 

1.00 

3 

0.00 

O.CO 

0 

-0.00000 

0.00000 

11.67 

2.79 

2.00 

3 

0.04 

0.02 

0 

-0.00000 

o.coooo 

11.67 

2.79 

3.00 

4 

-0.02 

-0.01 

0 

O.OOCOO 

o.coooo 

11.67 

2.79 

3.60 

4 

-0.06 

-0.04 

0 

-e.oooco 

o.ooooc 

11.67 

2.79 

0.0 

1 

0.0 

0.0 

0 

-0.00000 

0.00000 

11.99 

0.57 

1.00 

2 

-0.02 

-O.CO 

0 

-0.00000 

-o.coooo 

11.99 

0.57 

2.00 

3 

0.01 

0.00 

0 

-0.00000 

-o.coooc 

11.99 

0.57 

3.00 

V 

-0.00 

-0.00 

0 

-0.00000 

o.coooo 

11.99 

0.57 

3.30 

4 

-0.01 

-o.co 

0 

0.00000 

c.coooo 

11.99 

0.57 

0.0 

1 

0.0 

0.0 

0 

O.OOOUO 

o.coooo 

10.53 

-5.75 

0.91 

2 

-0.00 

0.00 

0 

0.00000 

o.coooo 

10.53 

-5.75 

0.0 

1 

0.0 

0.0 

0 

-O.OOOUO 

O.OOOCO 

11.30 

-4.04 

1.00 

2 

-e.oc 

0.00 

0 

0.00000 

0.00000 

11.30 

-4.04 

2.00 

3 

0.00 

-0.00 

0 

-0,000  00 

0.00000 

11.30 

-4.04 

2.03 

3 

0.00 

-0.00 

0 

0.00000 

0.00000 

11.30 

-4.04 

l 


H-7 


S2R-79-200 

Table  II  f age  -8- 

(p-Matrix) 


SOLUTIONS  OF  LOR AN  ECLAT IOFS  FOR  HORIZONTAL  POSITICN  AKO  VELOCITY  COORDINATES 


I 

K 

S  (KYOSt 

ITER 

ATE  (FT  ) 

CTE  (FT) 

1 

1 

0.0 

1 

0.0 

0.0 

L 

t 

l.LO 

3 

0.05 

-0.01 

t 

I 

2.20 

S 

0.03 

-0.00 

1 

1 

3.30 

6 

-0.11 

0.01 

t 

1 

4. AO 

7 

0.28 

-0.03 

1 

l 

5.50 

9 

0.21 

-0.02 

1 

1 

6.60 

11 

0.20 

-0.02 

l 

I 

7.70 

13 

0.22 

-0.02 

l 

1 

8.80 

15 

0.28 

-0.03 

1 

1 

9.90 

17 

0.37 

-0.03 

1 

l 

11.00 

19 

0.45 

-0.04 

t 

1 

11.19 

20 

-0.37 

0.03 

2 

1 

0.0 

1 

0.0 

0.0 

2 

1 

l.CO 

3 

0.00 

0.00 

2 

1 

2.00 

4 

-0.01 

-0.01 

2 

1 

3.00 

5 

0.02 

0.01 

2 

I 

A. 00 

6 

-0.03 

-0.02 

2 

1 

5.00 

7 

0.04 

0.03 

2 

1 

5. 42 

7 

o.oe 

0.05 

3 

1 

0.0 

1 

-o.oc 

-0.00 

3 

1 

0.66 

2 

-0.02 

-0.02 

1 

0.0 

1 

-0.00 

-0.00 

i 

1.00 

2 

0.00 

0.00 

l 

1.56 

2 

0.01 

0.02 

i 

0.0 

1 

-0.00 

-0.00 

t 

1.00 

3 

o.oc 

0.00 

l 

2.00 

3 

0.03 

0.04 

i 

3.00 

4 

-0.02 

-0.03 

i 

4.00 

5 

0.02 

0.02 

i 

4.29 

5 

0.03 

0.04 

i 

0.0 

1 

-O.OC 

-0.00 

l 

1.00 

3 

O.OC 

0.00 

l 

2.00 

3 

0.04 

0.02 

i 

3.00 

4 

-0.02 

-0.01 

i 

3.60 

4 

-0.06 

-0.04 

l 

0.0 

1 

0.0 

0.0 

l 

1.00 

2 

-0.02 

-0.00 

i 

2.09 

3 

0.01 

0.00 

i 

3.00 

4 

-o.oo 

-0.00 

i 

3.30 

4 

-0.01 

-9.00 

i 

0.0 

1 

0.0 

0.0 

t 

0.91 

2 

-0.00 

0.00 

t 

0.0 

1 

0.0 

0.0 

l 

1.00 

2 

-0.00 

0.00 

l 

2.00 

3 

O.CO 

-0.00 

l 

2.03 

3 

0.00 

-0.00 

ITER 

ATE  (MPh) 

CTE  IHPH) 

ROOT  IHPH) 

YOOT  (HPH) 

0 

o.oooco 

-0.00000 

7.95 

-8.98 

0 

-0.85570 

0.07196 

7.95 

-8.98 

0 

-1.68949 

0.13937 

7.95 

-8.98 

0 

-2  .49070 

0.2007C 

7.95 

-8.98 

0 

-3.25013 

0.25459 

7.95 

-8.58 

0 

-3.96048 

0.30013 

7.95 

-8.98 

0 

-4.61720 

0.33654 

7.95 

-8.98 

0 

-5.21656 

0.36336 

7.95 

-8.98 

0 

-5.75746 

0.38342 

7.95 

-8.58 

0 

-6.240C7 

0.38774 

7.95 

-8.58 

0 

-6.66574 

0.38552 

7.95 

-8.98 

0 

-6.73397 

0.38428 

7.95 

-8.58 

0 

O.OOOCO 

-0.00000 

11.53 

3.32 

0 

-O.388C0 

-0.20716 

11.53 

3.32 

0 

-0.78824 

-0.4  3418 

11.53 

3.32 

0 

-1.20022 

-0.68107 

11.53 

3.32 

0 

-1.623  25 

-0.94732 

11.53 

3.32 

0 

-2.05627 

-1.23193 

11.53 

3.32 

0 

-2.23974 

-1.35595 

11.53 

3.32 

0 

-O.OOOCO 

0  .00000 

10  .64 

5.56 

0 

-0.19319 

-0.21767 

10.64 

5.56 

0 

O.OOOCO 

-0.00000 

5.27 

10.78 

0 

0.09326 

0.06097 

5.27 

10.78 

0 

0.15061 

0.09031 

5.27 

10.78 

0 

O.OOOCO 

0.00000 

9.63 

7.16 

0 

-0.09101 

-0.25746 

9.63 

7.16 

0 

-0.177C4 

-0.54267 

9.63 

7.16 

0 

-0.25776 

-0.85530 

9.63 

7.16 

0 

-0.333C5 

-1.19420 

9.63 

7.16 

0 

-0.35344 

-1.29532 

9.63 

7.16 

0 

-O.OOOCO 

-0  .00000 

11.67 

2.79 

0 

-0.31214 

-0.27939 

11.67 

2.79 

0 

-0.61304 

-0.54644 

11.67 

2.79 

0 

-0.90026 

-0.79762 

11.67 

2.79 

0  • 

-1.06412 

-0.93850 

11.67 

2.79 

0 

O.OOOCO 

-0.00000 

11.99 

0.57 

0  . 

-0.27329 

-0.09572 

11.99 

0.57 

0 

-0.52859 

-0.16935 

11.99 

0.57 

0 

-0.766  51 

-0.22267 

11.99 

0.57 

0: 

-0.83388 

-0.23482 

11.99 

0.57 

0 

O.OOOCO 

-0.00000 

10.53 

-5.75 

0  . 

-0.03665 

0.11823 

10.53 

-5.75 

0 

0.00000 

-0.00000 

11.30 

-4.04 

0 

-0.09741 

0.13381 

11.30 

-4.04 

0 

-0.187S4 

0.27544 

11.30 

-4.C4 

0 

-0.190” 

0.27923 

11.30 

-4.04 
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TO:  C.  R.  Edwards 

FROM :  V .  S  chwab 

SUBJECT:  Computation  of  GDOP  for  St.  Mary's  River  Data  File 

The  formulas  used  to  compute  GDOP  for  the  St.  Mary's  River 
Data  File  are  summarized  below. 

1.  General  Formulas 

The  scalar  variable  GDOP  (geometric  dilution  of  £recision) 
depends  on  the  G-matrix  and  is  given  by 

GDOP  =  Itr(GGT)]1/2  (1) 

T 

in  which  G  denotes  the  transpose  of  the  matrix  G  and  tr(M)  denotes 
the  sum  of  the  diagonal  elements  of  M.  The  matrix  M=GGT  is  n  x  n 
where  n=2  or  3  according  as  2  or  3  LORAN  TDs  are  used  in  the  navi¬ 
gation  fix. 

In  the  case  of  2  TDs  the  G-matrix  is  given  by 

G  =  A-1  (2) 

where  A  is  2  x  2  and  of  rank  2. 

In  the  case  of  3  TDs  minimum  variance  estimation  is  used 
and  G  is  given  by 

T  -1  -1  t  -l 

G  *  (A  R  A)  A  R  .  (3) 

In  this  case  A  is  3  x  3  and  of  rank  3  and  the  3x3  matrix  R  is  the 
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covariance  matrix  of  the  noise  in  the  TD  measurements.  In  both  (2) 
and  (3)  the  elements  of  A  depend  only  on  the  relative  geometry  of 
the  navigat .  ,  master  station  and  slave  stations  and  the  nominal 
propagation  velocity  of  the  LORAN  signals. 

We  define  the  following  coordinates  which  enter  into  the 
definitions  of  the  elements  of  A  and  R: 


x  ,  y  horizontal  cartesian  position  coordinates 

of  the  master  station 

x.,y.  horizontal  cartesian  position  coordinates 

1  of  ith  slave  station 

x  ,y  horizontal  cartesian  position  coordinates 

p  p  of  a  known,  arbitrary  point 

v  nominal  propagation  velocity  of  the  LORAN 

signals 

R^  slant  range  from  master  station  to  arbi¬ 

trary  point 

R_.  slant  range  from  ifc^  slave  station  to  arbi¬ 

trary  point 

a. .  element  of  A-matrix 


The  formulas  for  R^,  R^,  aii 

R1  *  +  <VYm)2 

R2i  ’  [(V*i>2  +  <yp-yi> 


and  a 


i2 


1/2 

] 


are 


i-1.2,3 


(4) 

(5) 
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a.,  =  [ (x  -x.)/R,.  -  (x  -x  )/R,]/v,  i=l,2,3 
ll  p  i  2i  p  m  i 


ai2  "  I(Vyi>/R2i  "  <Vy»,/Rl,/V'  i=1'2'3 


2.  Formulas  for  GDOP  in  the  Case  of  Two  TPs 

In  the  case  of  two  TDs  eqs  (1)  and  (2)  give 


2  2  2  2  1/2 

GD0P2  =  ^all+a12+a21+a22^  /^alla22-a12a21^ 


3 .  Formulas  for  GDOP  in  the  Case  of  Three  TDs 

The  elements  of  the  noise  covariance  matrix,  R,  are 


rii  “  <R1  +  R2i>/Rl  '  i=1'2'3 


r^j  =  1.0,  l^j. 


The  determinant  of  R  is 


A  =  rllr22r33+  2r12r23r13"  r13r22_  r12r33  “  rllr23‘ 


The  elements  of  R  are 


(r22r33“r23)/A 


(r13r23_r33r12)/A 


(r12r23"r22r13)/A 
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r22  =  (rllr33“r13)/A 
r23  =  (r12r13“rllr23)/A 
r33  =  (rllr22_r12)/A 


-1 

-1 

r21 

“  r12 

-1 

-1 

r31 

"  r13 

-1 

-1 

r32 

“  r23 

The 

elements 

of  B  =  A1 

1  . 
R 

-1 

are 

-1 

-1 

-1 

bll 

“  allrll 

+ 

a21r21 

+ 

a31r31 

-1 

-1 

-1 

b12 

*  allr12 

+ 

a21r22 

+ 

a31r32 

-1 

-1 

-1 

b13 

"  allr13 

+ 

a21r23 

+ 

a31r33 

-1 

-1 

-1 

b21 

=  a12rll 

+ 

a22r21 

+ 

a32r31 

-1 

-1 

-1 

b22 

=  a12r12 

+ 

a22r22 

+ 

a32r32 

-1 

-1 

-1 

b23 

=  a12r13 

+ 

a22r23 

+ 

a32r33 

The  elements  of 

C  =  (A1 

i  _ 

R 

■1.. 

A)  are 

C11 

"  bllall 

+ 

b12a21 

+ 

b13a31 
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c12  "  blla12  +  b12a22  +  b13a32 

c2l  “  b21all  +  b22a21  +  b23a21 

c22  “  b21a12  +  b22a22  +  b23a32 

—1  T  —  1  —  1 

The  elements  of  C  =  (A  R  A) 


c22cll 

C22/A 

-cl2/t 

c2i/A 

°U/A 


C12C21 


The  elements  of  G  •  (A^R  A)  A  R  are 


gll  - 

Cllbll 

+ 

C12b21 

®12  = 

-1, 

cllb12 

+ 

c12b22 

*13  = 

cllb13 

+ 

c12b23 

921  “ 

c21bll 

+ 

c22b21 

1-5 


THE  JOHNS  HOPKINS  UNIVERSITY 

applied  physics  LABORATORY 

l*u*EL  MAftnAhO 


S2R-80-019 
Page  -6- 


922  =  C21b12  +  C22b22 
g23  =  C21b13  +  C22b23 
Finally 

1/2 

gdop3  =  [g^  +  g*2  +  gL23  +  g221  +  g22  +  (g> 

4.  The  Units  in  Which  GDOP  is  Expressed 

The  units  in  which  GDOP2  given  by  (8)  and  GDOP3  given  by 
(9)  are  expressed  depend  on  the  units  used  to  express  the  propa¬ 
gation  velocity,  v.  For  the  St.  Mary's  River  Data  File  the  position 
coordinates  are  expressed  in  kiloyards  and  propagation  velocity  in 
kiloyards/nanosec.  GDOP2  and  GDOP3  as  given  by  (8)  and  (9)  are  also 
then  in  kiloyards/nanosec.  However,  in  the  Data  File  it  was  desired 
to  express  GDOP  in  yds/nanosec,  so  the  numerical  values  given  by  (8) 
and  (9)  are  multiplied  by  1000  to  make  the  conversion. 


V.  Schwab 
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